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METHOD AND SYSTEM FOR TRACKING 
MULTIPLE REGIONAL OBJECTS BY 
MULTI-DIMENSIONAL RELAXATION 

RELATED APPLICATIONS 

This application is a continuation-in-part of U.S. patent 
application Ser. No. OS/404,024, filed Mar. 14, 1995 and 
issued Jul. 16, 1996 as U.S. Pat. No. 5,537,119, which is a 
continuation-in-part of U.S. patent application Ser. No. 
08/171327 filed Dec. 21, 1993, now U.S. Pat. No. 5,406, 
289. 

FIELD OF THE INVENTION 

The invention relates generally to computerized tech- 
niques for processing data obtained from radiation reflec- 
tions used to track multiple discrete object. 

BACKGROUND OF THE INVENTION 

The invention relates generally to computerized tech- 
niques for processing data obtained from radar to track 
multiple discrete objects. 

There are many situations where the courses of multiple 
objects in a region must be tracked. Typically, radar is used 
to scan the region and generate discrete images or "snap- 
shots" based on sets of returns or observations. In some 
types of tracking systems, all the returns from any one object 
are represented in an image as a single point unrelated to the 
shape or size of the objects. 'Tracking" is the process of 
identifying a sequence of points from a respective sequence 
of the images that represents the motion of an object. TTie 
tracking problem is difficult when there are multiple closely 
spaced objects because the objects can change their speed 
and direction rapidly and move into and out of the line of 
sight for other objects. The problem is exacerbated because 
each set of returns may result from noise as well as echoes 
from the actual objects. The returns resulting from the noise 
are also called false positives. Likewise, the radar will not 
detect all echoes from the actual objects and this phenomena 
is called a false negative or "missed detect* 1 error. For 
tracking airborne objects, a large distance between the radar 
and the objects diminishes the signal to noise ratio so the 
number of false positives and false negatives can be high. 
For robotic applications, the power of the radar is low and 
as a result, the signal to noise ratio can also be low and the 
number of false positives and false negatives high. 

In view of the proximity of the objects to one another, 
varied motion of the objects and false positives and false 
negatives, multiple sequential images should be analyzed 
collectively to obtain enough information to properly assign 
the points to the proper tracks. Naturally, the larger the 
number of images that are analyzed, the greater the amount 
of information that must be processed. 

While identifying the track of an object, a kinematic 
model describing the object's location, velocity and accel- 
eration may be generated. Such a model provides the means 
by which the object's future motion can be predicted. Based 
upon such a prediction, appropriate action may be initiated. 
For example, in a military application there is a need to track 
multiple enemy aircraft or missiles in a region to predict 
their objective, plan responses and intercept them. 
Alternatively, in a commercial air traffic control application 
there is a need to track multiple commercial aircraft around 
an airport to predict their future courses and avoid collision. 
Further, in these and other applications, such as robotic 
applications, may use radar, sonar, infrared or other object 
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detecting radiation bandwidths for tracking objects. In 
particular, in robotic applications reflected radiation can be 
used to track a single object which moves relative to the 
robot (or vice versa) so the robot can work on the object. 

Consider the very simple example of two objects being 
tracked and no false positives or false negatives. The radar, 
after scanning at time t A , reports objects at two locations in 
a first observation set. That is, it returns a set of two 
observations {o u , o 12 }. At time U it returns a similar set of 

10 two observations {o 21 , o 22 } from a second observation set. 
Suppose from prior processing that track data for two tracks 
T t and T 2 includes the locations at t 0 of two objects. Track 
T, may be extended through the points in the two sets of 
observations in any -of four ways, as may track T 2 . The 

15 possible extensions of Tj can be described as: {T t , o u , o 21 }, 
{T 1( o n , o 22 }, {T x , o 12 , o 21 } and {T,, o 12 , o 22 }. Tracks can 
likewise be extended from T 2 in four possible ways 
including, <[T 2 , o 12 , o 2i }. FIG. 1 illustrates these five (out of 
eight) possible tracks (to simplify the problem for purposes 

20 of explanation) . The five track extensions are labeled h u , 
h 12 , h 13 , h J4 , and h 21 wherein h u is derived from {Ti, o u , 
°2i K n i2 i s derived from {T t , o u , o 22 }, h 13 is derived from 
{T t , o i2 , o 21 }, h 14 is derived from {T A , o i2 , o 22 }, and h 21 is 
derived from {T 2 , o n , o 2A ). The problem of determining 

25 which such track extensions are the most likely or optimal 
is hereinafter known as the assignment problem. 

It is known from prior art to determine a figure of merit 
or cost for assigning each of the points in the images to a 
track. The figure of merit or cost is based on the likelihood 

30 that the point is actually part of the track. For example, die 
figure of merit or cost may be based on the distance from the 
point to an extrapolation of the track. FIG. 1 illustrates costs 
6 21 5 22 21 modeled target characteristics. The function to 
calculate the cost will normally incorporate detailed char- 

35 acteristics of the sensor, such as probability of measurement 
error, and lirack characteristics, such as likelihood of track 
maneuver. 

FIG. 2 illustrates a two by two by two matrix, c, that 

^ contains the costs for each point in relation to each possible 
track. The cost matrix is indexed along one axis by trie track 
number, along another axis by the image number and along 
the third axis by a point number. Thus, each position in the 
cost matrix lists the cost for a unique combination of points 

45 ^and a track, one point from each image. FIG. 2 also 
illustrates a {0, 1} assignment matrix, z, which is defined 
with the same dimensions as the cost matrix. Setting a 
position in the assignment matrix to "one** means that the 
equivalent position in the cost matrix is selected into the 

50 solution. The illustrated solution matrix selects the {h 14 , 
h 21 } solution previously described. Note that for the above 
example of two tracks and two snapshots, the resulting cost 
and assignment matrices are three dimensional. As used in 
this patent application, the term "dimension" means the 

55 number of axes in the cost or assignment matrix while size 
refers to the number of elements along a typical axis. The 
costs and assignments have been grouped in matrices to 
facilitate computation. 
A solution to the assignment problem satisfies two 

60 constraints — first, the sum of the associated costs for assign- 
ing points ito a track extension is rninimized and, second, if 
no false positives or false negatives exist, then each point is 
assigned to one and only one track. 
When false positives exist, however, additional hypotheti- 

65 cal track extensions incorporating the false positives will be 
generated. Further note that the random locations of false 
positives will, in general, not fit well with true data and such 
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additional hypothetical track extensions will result in higher 
costs. Also note that when false negative errors exist then 
the size of the cost matrix must grow to include hypothetical 
track extensions formulated with "gaps" (i.e.. data omissions 
where there should be legitimate observation data) for the 
false negatives. Thus, the second criteria must be weakened 
to reflect false positives not being as signed and also to 
permit the gap filler to be multiply assigned With hypo- 
thetical cost calculated in this manner then the foregoing 
criteria for nutumization will tend to materialize the false 
negatives and avoid the false positives. 

For a 3-dimensional problem, as is illustrated in FIG. 1, 
but with N x (initial) tracks, N 2 observations in scan 1, N 3 
observations in scan 2, false positives and negatives 
assumed, the assignment problem can be formulated as: 
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where "c" is the cost and "z" is a point or observation 
assignment as in FIG. 2. 

The minimiz ation equation or equivalently objective 35 
function [1.0] (a) specifies the sum of the element by 
element product of the c and z matrices. The summation 
includes hypothesis representations Z tl ^ with observation 
number zero being the gap filler observation. Equation 
[1.01] (b) requires that each of the tracks T v . . . .T N be 40 
extended by one and only one hypothesis. Equation [1.0] (c) 
relates to each point or observation in the first observation 
set and requires that each such observation, except the gap 
filler, can only associate with one track but because of the 
"less than" condition it might not associate with any track. 45 
Equation [1.01] (d) is like [1.0] (c) except that it is appli- 
cable to the second observation set. Equation [1.0] (e) 
requires that elements of the solution matrix z be limited to 
the zero and one values. 

The only known method to solve Problem Formulation 50 
[1.0] exactly is a method called "Branch and Bound." This 
method provides a systematic ordering of the potential 
solutions so that solutions with a same partial solution are 
accessible via a branch of a tree describing all possible 
solutions whereby the cost of unexamined solutions on a 55 
branch are incrementally developed as the cost for other 
solutions on the branch are determined. When the develop- 
ing cost grows to exceed the previously known minimal cost 
(Le., the bound) then enumeration of the tree branch termi- 
nates. Evaluation continues with a new branch. If evaluation 60 
of the cost of a particular branch completes, then that branch 
has lower cost than the previous bound so the new cost 
replaces the old bound. When all possible branches are 
evaluated or eliminated then the branch that had resulted in 
the last used bound is the solution. If we assume that 65 
N^N^N^n and that branches typically evaluate to half 
there full length, men workload associated with "branch and 



bound** is proportional to (nHV^!) 2 . This workload is 
unsuited to real time evaluation. 

The Branch and Bound algorithm has been used in past 
research on the Traveling Salesman Problem. Messrs. Held 
and Karp showed that if the set of constraints was relaxed by 
a method of Lagrangian Multipliers (described in more 
detail below) then tight lower bounds could be developed in 
advance of enumerating any branch of the potential solution. 
By initiating the branch and bound algorithm with such a 
tight lower bound, significant performance improvements 
result in thait branches will typically evaluate to less than half 
their full length. 

Messrs. Frieze and Yadagar in dealing with a problem 
related to scheduling combinations of resources, as in job. 
worker and work site, showed that Problem Formulation 
[1.0] applied. They further described a solution method 
based upon an extension of the Lagrangian Relaxation 
previously mentioned. The two critical extensions provided 
by Messrs. Frieze and Yadagar were: (1) an iterative proce- 
dure that permitted the lower bound on the solution to be 
improved (by "hill climbing" described below) and (2) the 
recognition that when the lower bound of the relaxed 
problem was maximized, then there existed a method to 
recover the solution of the non-relaxed problem in most 
cases using parameters determined at the maximum The 
procedures attributed to Messrs. Frieze and Yadagar are only 
applicable to the 3-dimensional problem posed by Problem 
Formulation [1.0] and where the cost matrix is fully popu- 
lated. However, tracking multiple airborne objects usually 
requires solution of a much higher dimensional problem. 

FIGS. 1 -and 2 illustrate an example where "look ahead" 
data from the second image improved the assignment accu- 
racy for the first image. Without the look ahead, and based 
only upon a simple nearest neighbor approach, the assign- 
ments in the first set would have been reversed. Problem 
Formulation (1.0] and the prior art only permit looking 
ahead one image. In the prior art it was known that the 
accuracy of assignments will improve if the process looks 
further ahead, however no practical method to optimally 
incorporate look ahead data existed. Many real radar track- 
ing problems involve hundreds of tracks, thousands of 
observations per observation set and matrices with dimen- 
sions in the range of 3 to 25 including many images of look 
ahead. 

It was also known that the data assignment problem may 
be simplified (without reducing the dimension of the assign- 
ment problem) by eliminating from consideration for each 
track those points which, after considering estimated limits 
of speed and turning ability of the objects, could not 
physically !>e part of the track. One such technique, denoted 
hereinafter the "cone method." defines a cone as a continu- 
ation of each previously determined track with the apex of 
the cone at the end of the previously defined track. The 
length of the cone is based on the estimated rnaximum speed 
of the object and the size of the arc of the cone is based on 
the estimate maximum turning ability of the object Thus, 
the cone defines a region outside of which no point could 
physically be part of the respective track. For any such 
points outside of the cones, an infinite number could be put 
in the cost matrix and a zero could be preassigned in the 
assignment matrix. It was known for the tracking problem 
that these elements will be very common in the cost and 
selection matrices (so these matrices are "sparse"). 

It was also known in the prior art that one or more tracks 
which are substantially separated geographically from the 
other tracks can be separated also in the assignment prob- 
lem. This its done by examining the distances from each 
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point to the various possible tracks. If the distances from one observations in a particular scan to which the constraint 

set of points are reasonably short only in relation to one set is related. In general, there is a constraint for each 

track, then they are assigned to that track and not further observation of the scan, wherein the constraint indi- 

considered with the remainder of the points. Similarly, if a cates the number of track extensions to which the 

larger group of points can only be assigned to a few tracks, 5 observation may belong. 

then the group is considered a different assignment problem. In particular, the first optimization problem is formulated. 

Because the complexity of assignment problems increases generated or defined as an M-dimensional assignment prob- 

dramatically with the number of possible tracks and the total lem wherein there are M constraint sets in the first collection 

number of points in each matrix, this partitioning of the of constraint sets (i.e., there are M scans being examined) 

group of points into a separate assignment problem and 10 and the first objective function minimizes a total cost for 

removal of these points from the matrices for the remaining assigning observations to various track extensions wherein 

points, substantially reduces the complexity of the overall terms are included in the cost, such that the terms have the 

assignment problem. figures of merit or costs for hypothesized combinations of 

A previously known Multiple Hypothesis Testing (MHT) assignments of the observations to the tracks. Subsequently, 
algorithm (see Blackman, Multiple-Target Tracking with 15 the formulated M-dimensional assignment problem is 
Radar Applications. Chapter 10, Artech House Norwood solved by reducing the complexity of the problem by 
MA, 1986) related to formulation and scoring. The MHT generating one or more optimization problems each having 
procedure describes how to formulate the sparse set of all a lower dimension and then solving each lower dimension 
reasonable extension hypothesis (for FIG. 1 the set \h n . . . optimization problem. That is, the M-dimensional assign- 
h 24 }) and how to calculate a cost of the hypothesis {T r o ly , 20 ment problem is solved by solving a plurality of optimiza- 
o 2k \ based upon the previously calculated cost for hypoth- tion problems each having a lower number of constraint sets, 
esis {T ( , o l7 }. The experience with the MHT algorithm, The reduction of the M-dimensional assignment problem 
known in the prior art, is the basis for the assertion that look to a lower dimensioned problems is accomplished by relax- 
ahead through k sets of observations results in improved ing the constraints on the points of one or more scans 
assignment of observations from the first set to the track. 25 thereby permitting these points to be assigned to more than 

In theory, the MHT procedure uses the extendable costing one track extension. In relaxing the constraints, terms having 

procedure to defer assignment decision until the accumu- penalty factors are added into the objective function thereby 

lated evidence supporting the assignment becomes over- increasing the total cost of an assignment when one or more 

whelming. When it makes the assignment decision it then points are assigned to more than one track. Thus, the 

eliminates all potential assignments invalidated by the deci- 30 reduction In complexity by this relaxation process is itera- 

sion in a process called "pruning the tree." In practice, the tively repeated until a sufficiently low dimension is attained 

MHT hypothesis tree is limited to a fixed number of gen- such that the lower dimensional problem may be solved 

erations and the overwhelming evidence rule is replaced by directly by known techniques. 

a most likely and feasible rule. This rule considers each track In one embodiment of the invention, each k-dimensional 

independently of others and is therefore a local decision rule. 35 assignment problem 2<k=iM, is iteratively reduced to a k- 1 

A general object of the present invention is to provide an dimensional problem until a 2-dimensional problem is speci- 

efficient and accurate process for assigning each point object fied or formulated. Subsequently, the 2-dimensional prob- 

in a region from multiple images to a proper track and then lem formulated is solved directly and a "recovery** technique 

taking an action based upon the assignments. is used to iteratively recover an optimal or near-optimal 

A more specific object of the present invention is to 40 solution to each k-dimensional problem from a derived 

provide a technique of the foregoing type which determines (k-1) dimensional problem k-2, 3, 4, . . .M. 

the solution of a k-dimensional assignment problem where In performing each recovery step (of obtaining a solution 

"k" is greater than or equal to three. to a k-dimensional problem using a solution to a (k-1)- 

„ tTiniAnv terror; txtt rcvrTT/'YNT' dimensional problem) an auxiliary function, is utilized. In 

SUMMARY OF THE INVENTION ^ ^ J neaf . optimal golution t0 

The present invention relates to a method and apparatus a ^dimensional problem, an auxiliary function, »F W , k=4, 

for tracking objects. In particular, the present invention 5 M IS specified and a region or domain is determined 

tracks movement or trajectories of objects by analyzing wherein this function is maximized, whereby values of the 

radiation reflected from the objects, the invention being reg ion determine the penalty factors of the (k-1)- 

especially useful for real-time tracking in noisy environ- 50 dimensional problem such that another 2-dimensional prob- 

ments. lem can be formulated which determines a solution to the 

In providing such a tracking capability, a region contain- k-dimensional problem using the penalty factors of the 

ing the objects is repeatedly scanned to generate a multi- (k-l)-dimensionai problem. 

plicity of sequential images or data observation sets of the Each Auxiliary function of both lower dimensional 

region. One or more points (or equivalently observations), in 55 problem penalty factors and a solution at the dimension k at 

each of the images or observation sets are detected wherein which the penalized cost function is solved directly 

each such observation either corresponds to an actual loca- (typically a 2-dimensional problem). Further, in 

tion of an object or is an erroneous data point due to noise. determining, for auxiliary function *¥ k determined, and 

Subsequently, for each observation detected, figures of merit utilized to identify the peak region. Thus, gradients are used 

or costs are determined for assigning the observation to each ^ f or eacn 0 f the approximation functions *F 3 ¥ 4 determining 

of a plurality of previously determined tracks. * Afterwards, penalty factors penalty factors until tm-1 is used in deter- 

a first optimization problem is specified which includes: mining the penalty factors for the M-l dimensional prob- 

(a) a first objective function for relating the above men- lem. Subsequently, once the M-dimensional problem is 
tioned costs to potential track extensions through the solved (using a 2-dimensional problem to go from an (M-l) 
detected observations (or simply observations); and 65 dimensional solution to an M-dimensional solution), one or 

(b) a first collection of constraint sets wherein each more of the following actions are taken based on the track 
constraint set includes constraints to be satisfied by the assignments: sending a warning to aircraft or a ground or sea 
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facility, controlling air traffic, controlling anti-aircraft or Other features and benefits of the present invention will 

anti-missile equipment, taking evasive action, working on become apparent from the detailed description with the 

one of the objects. accompanying drawings contained hereinafter. 

According to one feature of this first embodiment of the 

present invention, the following steps are also performed 5 BRIEF DESCRIPTION OF THE FIGURES 

before the step of defining the auxiliary function. Aprelimi- nG 1 is a h of or ^ gets generated by a 

nary auxiliary function is defined for each of the lower scan of a ion ^ possMe tnxks within &e images or 

dimensional problems having a dimension equal or one ^ sets according to me prior ^ 

greater than the dimension at which the penalized cost w ^ ^ ., t ... c ^ 

Lction is solved directly. The preliminary auxiliary tunc- 10 A nG * 2 cos * and ^signment matrices for the 

tion is a function of lower order penalty values and a *** sets 1 accordtn S t0 the P nor arL 

solution at the dimension at which the penalized cost func- FIG. 3 is a block diagram of the present invention. 

tion was solved directly. In determining a gradient of the FIG. 4 is a flow chart of a process according to the prior 

preliminary auxiliary function, step in the direction of the art for solving a 3-dimensionai assignment problem. 

gradient to identify a peak region of the preliminary auxil- 15 pjQ, 5 i s a fl ow cnart of a process according to the present 

iary function and determine penalty factors at the peak invention for solving a k-dimensionai assignment problem 

region. Iteratively repeat the defining, gradient determining. where "k" is greater than or equal to 3. 

stepping and peak determining steps to define auxiliary nG 6 £$ a h of va[im$ functions used t0 explain ^ 

functions at successively higher dimensions until the auxil- pre sent invention 

iary function at 6-18 (k-1) dimension is determined. In an 2 o ,™ „ . ' / « , , j. r - 

alternative second embodiment of the present invention, , FIG. 7 is another block diagram of the present invention 

instead of reducing the dimentiality of the M-dimensional f * soivm S the k-dimenszonal assignment problem where 

assignment problem by a single dimension at a time, a k 15 **** or e ^ ual t0 3 * 

pluraHty of dimensions are relaxed simultaneously. This new » is a flowchart describing the procedure for solvmg 
strategy has the advantage that when the M-dimensional 25 a n-dimensional assignment problem according to the sec- 
problem is relaxed directly to a 2-dimensional assignment ond embodiment of the invention, 
problem* then all computations may be performed precisely DETAILED DESCRIPTION OF THE 
without utilizing an auxiliary function such as embodi- PREFERRED EMBODIMENTS 
ment More particularly, the second embodiment solves the 

first optimization problem (i.e.. the M-dimensional assign- 30 Referring now to the other figures in detail wherein like 

ment problem) by specifying (i.e., creating, generating, reference numerals indicate like elements throughout the 

formulating and/or defining) a second optimization problem several views, FIG. 3 illustrates a system generally desig- 

Hie second optimization problem includes a second objec- nated 100 for implementing the present invention. System 

tive function and a second collection of constraint sets 100 comprises, for example, a radar station 102 (note sonar, 

wherein: 35 microwave, infrared and other radiation bandwidths are also 

a) the second objective function is a combination of the contemplated) for scanning a region which may be, for 
first objective function and penalty factors or terms example, an aerial region (in aerial surveillance 
determined for incorporating M-m constraint sets of the applications) or a work region (in robotic applications) and 
first optimization problem into the second objective generating signals indicating locations of objects within the 
function; 40 region. The signals are input to a converter 104 which 

b) the constraint sets of the second collection include only converts the signals to data points or observations in which 
m of the constraint sets of the first collection of each ( or false Positive) « represented by a single 
constraints, point The output of the converter is input to and readable by 

wherein 2Sm5M-l. Note that, once the second optimiza- a computer 106. As described in more detail below, the 
tion problemhas been specified orformulatecL an optimal or 45 computer 106 assigns the points to respective tracks, and 
near-optimal solution is determined and that solution is used displays the tracks and extrapolations of the tracks on 
in specifying (i.e.. creating, generating, formulating and/or a monitor 110. Also, the computer 106 determines an 
defining) a third optimization problem of M-m dimensions appropriate action to take based on the tracks and track 
(or equivalent^ constraint sets). The third optimization extensions. For example, in a commercial application at an 
problem is subsequently solved by decomposing ii using the 50 airport, the computer can determine if two aircraft being 
same procedure of this second embodiment as was used to tracked are on a collision course and if so, signal a trans- 
decompose the first optimization problem above. Thus, a ***** 112 t0 wara each or tf a scheduled take-off 
plurality of instantiations of the third optimization problem wm P° se the risk of collision. <*elay the take-off. For a 
are specified, each successive instantiation having a lower military application on a ship or base, the computer can 
number of dimensions, until an instance of the third opti- 55 determine subsequent coordinates of enemy aircraft and 
mization problem is a two dimensional assignment problem send me coordinates to an antiaircraft gun or missile 120 via 
which can be solved directly. Subsequently, whenever an a communication channel 122. In a robotic application, the 
instance of the third optimization problem is solved, the computer controls the robot to work on the proper object or 
solution is used to recover a solution to the instance of the portion of the object 

first optimization problem from which this instance of the 60 The invention generates k-dimensional matrices where k 

third optimization was derived. Thus, an optimal or near- is the number of images or sets of observation data in the 

optimal solution to the original first optimization problem look ahead window plus one. Then, the invention formulates 

may be recovered through iteration of the above steps. a k-dimensional assignment problem as in [1.0] above. 

As mentioned above, the second embodiment of the The k-dimensional assignment problem is subsequently 

present invention is especially advantageous when m==2, 65 relaxed to a (k-l)-dimensional problem by incorporating 

since in this case all computations may be performed one set of constraints into the objective function using a 

precisely and without utilizing auxiliary functions. Lagrangian relaxation of this set Given a solution of the 
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(k-l)-dimensional problem, a feasible solution of the observation with k-1 gap fillers, e.g., Zq Ol o-. °/i p ^0. 
k-dimensional problem is then reconstructed. The Thus, the resulting generalization of Problem Formu- 
dimensional problem is solved in a similar manner, and the lation 1 1.0 1 without the "less than" complication within 
process is repeated until it reaches the two-dimensional the constraints is the following k-Dimensional Assign- 
problem that can be solved exactly. The ideas behind the 5 ment Problems in which k£3: 
Lagrangian relaxation scheme are outlined next. 

Consider me integer programming problem Minimize- V u ^ ^ 

Minimize v(z) =c r z u =o 1 * 1 * 

Subject to: Az ^b 10 
Bz ^d 

z, is an integer for i el ' ^ , i=0 
1 1.0.1 1 

where the partitioning of the constraints is natural in some n { s t 

sense. Given a multiplier vector u^O, the Lagrangian relax- 15 X — X 2 

ation of [1.0.1] relative to the constraints Bz^d is defined to v-i =0 ^+i =0 '* =0 
be 

<J><J>TT 

Subject to: Az^b and , = 2, ...*-i 

z t is an integer for i el 20 Wj Ni l 

|1:0.2| £ u^-U 

If the constraint set Bz^d is replaced by Bz=d, the nonne- 'i=o f*-i=o 

gativity constraint on u is removed. =C T z+u r (Bz-d) is the z, x ... i k e (0, l ), v / a n = l * 

Lagrangian relative to the constraints Bz^d, and hence the ^ 
name Lagrangian relaxation. The following fact gives the 

relationship between the objective functions of the original whcre c and z are similarly dimensioned matrices represent- 

and relaxed problems in S costs and hypothetical assignments. Note, in general, for 

FACT A.l. If i isan optimal solution to f 1.0.1], then tracing problems these matrices are sparse. 

<D£v(z) for ail u*0. If an optimal solution z of 11.0.2] is Ate formulating the k-dimensional assignment problem 

7 -i; x I maoT *u A * f"""* .1*-/ ii An 30 asm 1.1 , the present invention solves the resulting prob- 

feasible for [1.0.2], then z is an optimal solution for [1.0.1] , 1 ' * *u — ♦ • j w j - it« 

. J v lem so as to generate the outputs required by devices 110. 

and <PAigor ltnm 112 m and m Fof each observation set Cj xece ived from 

converter 104 at time t, where i=l . . . , Nj, the computer 

k k=o processes 0 ( in a batch together with the other observation 

35 sets 0^ +1 O t and the track T,.*, i.e., T, is the set of all 

- <^^m * tracks that have been defined up to but not including O,. 

converging to the solution u of Maximize <4>*0} and a (Note , bold type designations refer to the vector of elements 

corresponding sequence of feasible solutions {Z*} t =0of for me indicated ^ Le „ me set of ^ observations in the 

[1.0.1] as follows: scan 0f ex i stm g at ^ eta ) The result of this 

1. Generate in initial approximation u 0 . 40 processing is the new set of tracks T iWt+l and a set of cost 

2. Given u*, choose a search direction s* and a search weighted possible solutions indicating how the tracks might 
distance a k a k k *¥ k extend to the current time t,. At time t, +1 the batch process 

estimate u k by u Jt+i =u Jt +a it k is repeated using the newest observation set and deleting the 

3. Given u k+l and a feasible solution % M (u k+l ) of [1.0.2], oldest. Thus, there is a moving window of observation sets 
recover a feasible solution z k+l (u* +1 ) of the integer 45 which is shifted forward to always incorporate the most 
programming problem [1.0.1). recent observation set. The effect is that input observation 

4. Check the termination criteria. If the algorithm is not se * f u reused for ^} °y cles and ° n Nervation 
finished, set k=k+l and return to Step 2. Otherwise. f et . s k ' th re , u f e each *• ervaUon Wlthm the observation set 
terminate the algorithm. is migrated into a track. 

If z is an optimal solution of f 1.0. 1]. then 50 m ?' \ m f" tes various f* ocess ^ '"f ^ 

/-\ < / \ receipt of each observation set. Except for the addition of the 

n , k . ~ } Z . A . * r **^- » k-dimensional assignment solving process 300 and the 

Since the optimal solution 2=z(u) of [1.0.2] is usuaUy not a modification to scoring process 154 t0 build ^ struc mres 

feasible solution of [1.0.1) for any choice of the multipliers, suitable for process 300 ^ processes m mGt 7 « based 

55 upon prior art. The following processes 150 and 152 extend 

* previously defined tracks h^ A based on new observations. 

* ot a Gate formulation and output process 156 determines, for 

each of the previously defined tracks, a zone wherein the 
track may potentially extend based on limits of velocity, 

<p 60 maneuverability and radar precision. One such technique to 

a a accomplish this is the cone method described previously. 

The definition of the zone is passed to gating process 150. 
When a new observation set O, is received, the gating 

p =p=k in z*; i* in problem [1.1] below indicates that process 150 will match each member observation with the 

the track extension represented by z* x i* includes a 65 zone for each member of the hypothetical set h/.j. After all 

false positive in the p"* observation set. Note that this input observations from O, are processed the new hypothesis 

implies that a hypothesis be formed incorporating an set h, is generated by extending each track of the prior set of 
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hypothetical tracks h, , either with missed detect gap fillers observation identifiers, and passes them to k-dimensional 
or with' all new observation elements satisfying the track's assignment problem solving process 300. 
zone This is a many to many matching in that each The assignment solving process 300 is described below, 
hypothesis member can be extended to many new observa- Its output is simply the list of assignments which constitute 
tions and each new observation can be used to extend many 5 the most likely solution of the problem described by Equa- 
hypotheses. It however, is not a full matching in that any tion [1.1]. Note that both gate formulation and output 
hypothesis will neither be matched to all observations nor processes 156 and 162 use (at different times) the list or 
vice versa. It is this matching characteristic that leads to the assignments to generate the updated track history T, to 
sparse matrices involved in the tracking process. eliminate or prune alternative previous hypotheses that are 
Subsequently, gating 150 forwards the new hypothesis set h £ 10 prohibited by the actual assignments m the list, and subse- 
to filtering process 152. Filtering process 152 determines a quently to output any required data. Also note that when one 
smooth curve for each member of h r Such a smooth curve of the gate formulation and output processes 156 and 162 
is more likely than a sharp turn from each point straight to accomplish these tasks, the process will subsequently gen- 
the next point. Further, the filtering process 152 removes erate and forward the new set of gates for each remaining 
small errors that may occur in generating observations. Note is hypothesis and the processes will then be prepared to receive 
that in performing these tasks, the filtering process 152 the next set of observations. In one embodiment, the loop 
preferably utilizes a minimization of a least squares test of described here will generate zones for a delayed set of 
the points in a track hypothesis or a Kalman Filtering observations rather than the subsequent set This permits 
approach. processes 156 and 162 to operate on even observation sets 
As noted above, the foregoing track extension process 20 while the scoring step 154 and k-dimensional solving pro- 
requires knowledge of a previous track. For the initial cess operate on odd sets of observations, or vice versa, 
observations, the following gating process 158 and filtering The assignment solving process 300 permits the present 
process 160 determine the "previous tract* based on initial invention to operate with window size of k-1 for k=3. The 
observations In determining the initial tracks, the points upper limit on k depends only upon the computational power 
from the first observation set form the beeinning points of all 25 of the computer 106 and the response time constraints of 
possible tracks. After observation data from the next obser- system 100. The k-1 observation sets within the processing 
vation set is received, sets of simple two point straight line window plus the prior track history result m a k-dimensional 
tracks are defined. Then, promotion, gate formulation, and Assignment Problem as described by Problem Formulation 
output step 162 determines a zone in which future exten- [1.1]. The present invention solves this generalized problem 
sions are possible. Note that filtering step 160 uses curve 30 including the processes required to consider false positives 
fitting techniques to smooth the track extensions depending and negatives, and also the processes required to consider 
uponthe number of prior observations that have accumu- sparse matrix problem formulations, 
lated in each hypothesis. Further note that promotion, gate L A FIRST EMBODIMENT OF THE k-DIMENSIONAL 
formulation and output process 162 also determines when ASSIGNMENT SOLVER 300 

sufficient observations have accumulated to form a basis for 35 In describing a first embodiment of the k-dimensional 
promoting the track to processes 150 and 152 as described assignment solver 300. it is worthwhile to also discuss the 
3 bove process of FIG. 4 which is used by the solver 300. FIG. 4 
The output of each of the filtering processes 152 and 160 illustrates use of the Frieze and Yadagar process as shown in 
is a set of hypothetical track extensions. Each such extension prior art for transforming a 3-dimensional assignment prob- 
contains the hypothetical initial conditions (from the previ- 40 lem into a 2-dimensional assignment problem and then use 
ous track), the list of observations incorporated in the a hill climbing algorithm to solve the 3-dimensional assign- 
extension, and distance between each observation and the ment problem. The solver 300 uses a Lagrangian Relaxation 
filtered track curve. Scoring process 154 determines the technique (well known in the art) to reduce the dimension of 
figure of merit or cost of an observation being an extension an original k-dimensional assignment problem (k>3) down 
of a track. In one embodiment the cost is based on the 45 to a 3-dimensional problem and then use the process of FIG. 
above-mentioned distance although the particular formula 4 to solve the 3-dimensional problem. Further note that the 
for determining the cost is not critical to the present inven- Lagrangian Relaxation technique is also utilized by the 
tion. A preferred formula for determining the cost utilizes a process of FIG. 4 and that in using this technique the 
negative log likelihood function in which the cost is the requirement that each point is assigned to one and only one 
negative of the sum of: (a) the logs of the distances nor- 50 track is relaxed. Instead, an additional cost, which is equal 
malized by sensor standard deviation parameters, and (b) the to a respective Lagrangian Coefficient "u". is added to the 
log likelihoods for events related to: track initiation, track cost or objective function [1.0] (a) whenever a point is 
termination, track maneuver, false negatives and false posi- assigned to more than one track. This additional cost can be 
tives. Note that track maneuvers are detected by comparing picked to weight the significance of each constraint violation 
the previous track curve with the current extension. Further 55 differently, so this additional cost is represented as a vector 
note that some of the other events related to, for example, of coefficients u which are correlated with respective obser- 
false negatives and false positives are detected by analyzing vation points. Hill climbing wifi then develop a sequence of 

the relative relationship of gap fillers in the hypothesis. Lagrangian Coefficients sets designated (u 0 u y , . . . 

Thus, after determining that one of these events occurred, a u p ). That correspond to an optimum solution of the 
cost for it can be determined based upon suitable statistics 60 2-dimensional assignment problem. The assignments at mis 

tables and system input parameters. The negative log like- optimum solution are then used to "recover" the assignment 

lihood function is desirable because it permits effective solution of the 3-dimensional assignment problem, 

integration of the useful components. Copies of the set of In step 200 of FIG. 4, initial values are selected for the u 0 

hypothetical track extensions which are scored are subse- coefficients. Because the Lagrangian Relaxation process is 
quently passed directly to one of the gate formulation and 65 iterative, the initial values are not critical and are all initially 

output steps 156 and 162. Note that the scoring process 154 selected as zero. In step 202, these additional costs are 

also arranges the actual scores in a sparse matrix based upon applied to the objective function [1.0] (a). With the addition 



5,959.574 

13 14 

of the costs "u", the goal is still to assign the points which 

imnimize the total cost. This transforms Equation 1 1.0] (a), * * 

written for k=3 and altered to exclude mechanisms related to 

false positives and negatives, into objective function |2.1| 

(a) In the first iteration it is not necessary to consider the "u" 5 

matrix because all "u" values are set to zero. To relax the 

requirement that each point be assigned to one and only one 

track, the constraint Equation 1 1.0] (d) is deleted, thereby coefficients whose corresponding Problem Formulation 
permitting points from the last image to be assigned to more 10 ' 2 ' 2 ' P roblera solution is a cost value closer to * c *** of 
than one track. Note that while any axis can be chosen for 

relaxation, observation constraints are preferably relaxed. dimensional assignment problem requires solving the Equa- 
The effect of this relaxation is to create a new problem which tion 2.2 problem corresponding to u p . 
must have the same solution in the first two axes but which ^ st£p 206 ^ WQ dimensional problem is solved 
can have a differing solution in the third axis. The result is dkectly using a technique to those skUled m me art 

constraints [2.1] (b-d). &uc ^ as R everse Auction for the corresponding cost and 

^ [2 ^ solution values. This is the evaluation of one point on the 

(a) Minimize: £ £ {(c,^ - % * surface or for the fet iteration W o 

' l=0i2=0 20 Thus, after this first iteration, the points have been 

n 2 Af 3 assigned based on all "u" values being arbitrarily set to zero. 

(b) Subject to: £ £ z, |i2 , 3 = 1, i, = 1, Because ttite "u" values have been arbitrarily assigned, it is 

' 2=I * =I unlikely that these assignments are correct and it is likely 

Nl ^ 25 that further iterations are required to properly assign the 

( c) £ Y^zi&ts z i, h ~ i Ni points. Step 208 determines whether the points have been 

'i =l '3=» properly assigned after the first iteration by determining if 

, cjn u y, T for this set of assignments whether a different set of "u" 

values could result in a higher total cost. Thus, step 208 is 
30 implemented by determining the gradient of objective func- 
tion 1 2.2 1 (a) with respect to u y . If the gradient is substan- 
Step 204 then generates from the 3-dimensional problem tialiy non-zero (greater than a predetermined limit) then the 
described by Problem Formulation [1.0 J a new assignments are not at or near enough to the peak of the <t>u 
2-dimensional problem formulation which will have the 35 

same solution for the first two indices. As Problem Formu- ^ y+l . . „ , . t t 

lation 12.1] has no constraints on the 3„ axis, any value Hm M ^ Step 212 determines the u values based 
within a particular 3 rd axis can be used in a solution, but u P on the u j val " es < *? ^ecuon resulting from protecting 
using anything other than the minimum value from any 3-rd the P r ™ s S"*** f <> *f U doma ^ a 
axis has the effect of increasing solution cost. Conceptually, " ™* soluto « value of 2-dimensional prob em is the set 
the effect of step 204 is to change the 3-dimensionai arrays of coefficients that minimize the 2-dimensional problem and 
in Problem Formulation [2.11 into 2-dimensional arrays as the actual cost at the minimum. Those coefficients aug- 
shown in ProblemFormulation [2.2] and to generate the new mented b > r the coefficients stored in ny 2 permit the evalu- 
2-dimensional matrix m t * 2 defined as shown in Equation 45 adon (but not the minimization) of the cost term m Problem 
j2 3] 1 Formulation [2.1]. These two cost terms are lower and upper 

bounds on the actual imnimized cost of the 3-dimensional 
«i *2 ( Muc ) I 2 - 2 ! problem, and the difference between them in combination 

( a) Minimize: Y Y (c ilhli - u h )L, 2 with the gradient is used to compute the step size. 

1^o£o ( vij 3 1 so With this new set of "u w values, steps 202-210 are 

repeated as a second iteration. Steps 212 and 202-210 are 

(b) Subject to- V V = 1 i, = l N x repeated until the gradient as a function of u determined in 

J 3 * ' 1 step 208 is less than the predetermined limit. This indicates 

that the up values which locate the peak area of the 0 u 

V y zi i ^ l, k = 1 N 2 55 surface are determined and that the corresponding Problem 

ijA & 1 3 Formulation [2.2] has been solved. Step 214 will attempt to 

use the assignments that resulted from this particular 
2-dimensional assignment problera to recover the solution of 

= Mia: arg minimize {c^u - u k jv*» * 2 60 the 3-dimensional assignment problem as described below. 

If the limiit was chosen properly so that the "u" values are 
close enough to the peak, this recovery will yield the proper 
set of assignments that rigidly satisfies the constraint that 
The cost or objective function for the reduced problem as each poinlt be assigned to one and only one track. However, 
defined by [2.2] (a), if evaluated at all possible values of u 65 if the V values are not close enough to the peak, then the 
is a surface over the domain of U 3 . This surface is referred limit value for decision 210 is reduced and the repetition of 
to as <I>u steps 212 and 202-210 is continued. 



to 

(d) ^€{0,1} VZi li2 



and j=Z...L-l 
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Step 214 recovers the 3-dimensional assignment solution 

by using the assignment values determined on the last ^ ^ [3.1] 

iteration through step 208. Consider the 2-dimensional z < a > Minimis: Vt tf)= ^-X^-'X.-'* 

assignment matrix to have l*s in me locations specified by * 1 ' '** 

the list LMa,. b s f^. If the 3-dimensional z matrix is 5 ^ f § ,> u „ . u ... Nl 

specified by placing l's at the location indicated by the list fa fa 1 ' * 

LyKEf, b,. rry,)^ then the result is a solution of Problem ^ ^ ^ 

Formulation [2.1]. Let I^m^Atei be me fanned by (c) E "* Z Z -Z^i- 1 * = 1 

the third index. If each member of L 3 is unique then the L 2 io 'I s0 v+j-* 

solution satisfies the third constraint so it is a solution to for ^ = x Nj 
Problem Formulation [1.0]. When this is not the case, 
recovery determines the minimal substitutions required 

within list L 3 so that it plus L A will be a feasible solution. i.e., ^ y y N k _ x ^ ~\ h = i N k 

a solution which satisfies the constraints of a problem 15 "i k „ x l '"' k * 

formulation, but which may not optimize the objective . . 

^, (e) & t €{0, If, for all h 

function of the problem formulation. This stage of the l "' k 
recovery process is formulated as a 2-dimensional assign- 
ment problem: Form a new cost matrix [c v ]" Ii/=1 where 20 and where c* is the cost matrix [c* . . .J which is a function 
for H • • • N 4 and the N, term is the total number of the between the observed point z* and the 

of cost elements in the selected row of the 3-dimensional smoothed track determined by the filtering step, and \ k is the 
cost nmtrix. Attempt to solve this 2-dimensional problem for cost f unct i on . xhi s set Q f equations is similar to the set 
the minimum using two constraints sets. If a feasible solu- presented in Problem Formulation [1.1] except that it 
tion is found then the result will have the same form as list 25 indudes tiie subscri p t m(i superscript k notation. Thus, in 
L x . Replace the first set of indexes by the indicated (a,, bj so lving the M-dimensional Assignment Problem the inven- 
pairs taken from list L x and the result will be a feasible tion reduces ^ problem to an M-l dimensional Assign- 
solution of Problem Formulation [1.1]. If no feasible solu- ment problem and then to an N-2 dimensional Assignment 
tion to the new 2-dimensional problem exists then further 3Q problem, etc. Further, the symbol ke{3. . . . M] customizes 
effort to locate the peak of O u Problem Formulation [3.1] to a particular relaxation level 

That is. the notation is used to reference data from levels 

* relatively removed as in c* +1 are the cost coefficients which 

existed prior to this level of relaxed coefficients c*. Note that 
0 * 35 actual observations are numbered from 1 to N t where N ( is 

the number of observations in observation set i. Further note 
that the added 0 observation in each set of observations is the 
unconstrained "gap filler." This element serves as a filler in 
* 40 substituting for missed detects, and a sequence of these 

* elements including only one true observation represents the 

possibility that the observation is a false positive. Also note 
that by being unconstrained a gap filler may be used in as 
many hypotheses as required, 

i 45 

While direct solution to [3.1] would give the precise 
assignment, the solution of k-dimensional equations directly 

by the set of observations {o=li=l M-l}. is not valid. for lar S e k is t0 ° ""^ consuming for practice. 

Because the matrix is sparse the list of cost elements is ™ us - me P resent invention solves P roblem 
stored as a packed list and then for each dimension of the 50 The following is a short description of many aspects of the 
matrix, a vector of a variable length list of pointers to the present invention and includes some steps according to the 
cost elements is generated and stored. This organization P™* The step in solving the problem indirectly is 
means that for a particular observation o= the j* list in the i'* t0 reduce * e complexity of the problem by the previously 
vector will be a list of pointers to all hypotheses in which o 55 known ™ d ^cussed lagrangian Relaxation technique, 
y participates. This structure is further explained in the According to the Lagrangian Relaxation technique, the 
Mowing section dealing with problem partitioning. ab f lute r^ufrement that each point is assigned to one and 
& ° r only one track is relaxed such that for some one image. 

The objective of the assignment solving process is to points can be assigned to more than one track. However, a 
select from the set of all possible combinations of track ^ penalty bused on a respective Lagrangian Coefficient u* is 
extensions a subset that satisfies two criteria. First, each added to the cost function when a point in the image is 
point in the subset of combinations should be assigned to assigned to more than one track. The Lagrangian Relaxation 
one and only one track and therefore, included in one and technique reduces the complexity or "dimension" of the 
only one combination of the subset, and second, the total of formulation of the assignment problem because constraints 
the scoring sums for the combinations of the subset should 65 on one observation set are relaxed. Thus, the Lagrangian 
be niiniraized. This yields the following M-dimensional Relaxation is performed iteratively to repeatedly reduce the 
equations where k=M: dimension until a 2-dimensional penalized cost function 
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problem results as in Problem Formulation |2.1|. This 

2- dimensional problem is solved then directly by a previ- function * a a 
ously known technique such as Reverse Auction. The penal- 
ized cost function for the 2-dimensional problem defines a * 
valley or convex shaped surface which is a function of 5 
various sets of {u*lk=3, . . . ,M} penalty values and one set 
of assignments for the points in two dimensions. That is, for 

each particular u 3 there is a corresponding 2-dimensional «* 

penalized cost function problem and its solution. Note that 

, . . . i. j ^ £ 10 gradient of the * 
the solution of the 2-dimensionai penalized cost function 

problem identifies the set of assignments for the particular u 3 u* u 
values that minimize the penalized cost function. However, 
these assignments are not likely to be optimum for any 

higher dimensional problem because they were based on an 15 results in a new 3-dimensional problem and requires the 

initial arbitrary set of u* values. Therefore, the next step is 2-dimensional hill climbing based upon new 2-dimensional 

to determine the optimum assignments for the related problems. At the peak of the 3-dimensional H* to a complete 

3- dimensional penalized cost function problem There exists solution. This process is repeated iteratively until the u' 
a 2-dimensional hill shaped function cl> values that result in a corresponding solution at the peak of 

20 the highest order 4* <t> 



u k \k > 3[ 



values originally assigned. Then, the gradient of the hill- 
shaped 4> 



25 



30 



35 



previously selected for the one point on the hill 
(corresponding to the minimum of the penalized cost <i> 



U 3 



40 



which the corresponding problem will result in the peak of 
the 



45 



50 



55 



Coefficient u M penalty values initially equal to zero. The 
Lagrangian Coefficients associated with the M"* constraint 
set are a N M+l element vector. The reduction of this 
M -dimensional problem in step 322 to a M-l dimensional 
problem uses the two step process described above. First, a 
penalty based on the value of the respective u M coefficient 
is added to the cost function when a point is assigned to more 
than one track and then the resultant cost function is mini- 
mized. However, during the first iteration, the penalty is zero 
because all u M values are initially set to zero. Second, the 
requirement that no point from any image can be assigned to 
more than one track is relaxed for one of the images. In the 
characteristic. The next task is to recover the solution of the ^ extreme this would allow a point from the relaxed image to 
proper u 4 values for the 4-dimensional problem. The fore- associate with every track. However, the effect of the 
going hill climbing process will not work again because the previous penalty would probably mean that such an asso- 
foregoing hill climbing process when required to locate the ciation would not minimize the cost The effect of the two 
4-dimensional u 4 values for the peak of <J> 3 exact definition s t e ps in combination is to remove a hard constraint while 
of the O 3 case of <I> 2 O 3 65 adding the penalty to the cost function so that it operates like 

a soft constraint. For step 322 this two step process results 
the iteration can result in a greater than approximation of O 3 in the following penalized cost function problem with fc=M: 
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-continued 



i/iZ-s* 4 *-- 



Subject to: .. . £ j*" 1 , (jL = 1 ii = 1 Ni 



Z-ZZ-£«.-' 

<1=0 ^i=0^ +1 *0 * 



0) Subject to: V ... £ a*, Jjt = 1 n = 1. ... N t 



10 



15 



for i ; « 1, ... AT, 
and y = 2, ... X - 1 

jf,„. lt 6 10.1) for all ii. .I t - 



*— 4 i, -=o 



Z-ZZ-S^- 

for ij=U...Nj 
and y = 2,...Jt-l 

<.„ <4 €{0,1} for all i,.../ t 



25 



Because the constraint on single assignment of elements 
from the last image has been eliminated, a M-l dimensional 
problem can be developed by eliminating some of the 
possible assignments. As shown in Equations [3.3]* this is 
done by selecting the smallest cost clement from each of the 
M th axis vectors of the cost matrix. Reduction in this manner 
yields a new. lower order penalized cost function defined by 
Equations [3.3] which has the same minimum cost as does 
the objective function defined by [3.2] (a) above. 

The cost vectors are selected as follows. Define an index 



array m* £ . • • and new cost array c 



-iby: 



20 Assignment Problem [3.1] and Problem Formulation [3.4] 
differ only in the dimension M vs. M-l, respectively. An 
optimal solution to [3.4] is also an optimal solution to 
equation [3.2]. This relationship is the basis for an iterative 
sequence of reductions indicated by steps 320-332 through 
330-332 and 200-204 in which the penalized cost function 
problem is reduced to a two dimensional problem. As these 
formula will be used to describe the processing at all levels, 
the lowercase k is used except where specific reference to 
30 the top level is needed. In step 206. the 2-dimensional 
penalized cost function is solved directly by the prior art 
Reverse Auction technique. Each execution of 206 produces 
two results;, the set of z 2 values that inMrnize the problem 
and the cost that results v 2 when these z values are substi- 
tuted into the objective function [3.4] (a). 

In step 208, according to the prior art. solution z 2 values 
are substituted into the 2-dimensional derivative of the <J> 2 
U should be adjusted so as to perform the hill climbing 
40 function. As was previously described the objective is to 
produce a sequence of u 3 , values which ends when the U 3 p 
value in the domain of the peak of the 4> 



35 



wf lmii = Mini arg minimize^ - «fj i k = 0, 1 N k \ 



[33] 



45 



4:U 



= £mia{0,4... 0 -«* t ) 
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The resulting M-l dimensional problem is (where fc=M): 



[3.4] 



55 



to the peak of <J> 2 flow moves to step 214 which will attempt 
to recover the 3 -dimensional solution as previously 
described. When further adjustment is required then the flow 
progresses to step 212 and the new values of u* are com- 
puted. At the 2-dimensional level the method of the prior art 
could be u;sed for the hill climbing procedure. However, it is 
not practical to use this prior art hill climbing technique to 
determine the updated Lagrangian Coefficients u* or the 
Max on the next (or any) higher order <J> 
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it is not possible to use the prior art hill climbing to hill climb 
to the peak of & k definition of <3> t functions. As the solution 
of 4>\ in that higher order u values are not yet optimized, its 
value can be larger than the true peak of Ok lower bound on 
v A and it can not be used to recover the required solution. 



10 



values which result in z values that are closer to the proper 
z values for the highest k-dimension. This hill climbing 
process (which is different than that used in the prior art of 
FIG. 4 for recovering only the three dimensional solution) is 
used iteratively at all lower dimensions k of the 
M-dimensional problem (including the 3-dimensional level 
where it replaces prior art) even when k is much larger than 
three. FIG. 6 helps to explain this new hill climbing tech- 
nique and illustrates the original k-dimensional cost function 
Vjt of Problem Formulation |3.1|. However, the actual 
k-dimensional cost surface v* defined by [3.1 1 (a) comprises 
scalar values at each point described by k-dimensional 
vectors and as such can not be drawn. Nevertheless, for 
illustration purposes only, FIG. 6 ignores the reality of 
ordering vectors and illustrates a concave function v^z*) to 
represent Equation 3.1. The \ k surface is illustrated as being 
smooth to simplify the explanation although actually it can 
be imagined to be terraced. The goal of the assignment 
problem is to find the values of z*; these values minimize the 
k-dimensional cost function v^ 

For purposes of explanation, assume that in FIG. 6, k=4 
(the procedure is used for all k£3). This problem is reduced 
by two iterations of Lagrangian Relaxation to a 
2-dimensional penalized cost function <p <;: ~ 1)i . This cost 
function, and all other cost functions described below are 
also non-smooth and continuous but are illustrated in FIG. 
6 as smooth for explanation purposes. Solving the 
problem results in one set of z 2 assignments and the value of 
<&<*-i) at the point u 2 0 . A series of functions .... 
^<*-ft each generated from a different u 3 are shown. The 
particular series illustrates the process of locating the peak 
of <t>/*_ lv The 2-dimensional penalized cost functions 
^<*-i\ /. . , can be solved directly. Each such 

solution provides the information required to calculate the 
next u 3 value. Each iteration of the hill climbing improves 
the selection of u 3 values, i.e., yields 4> 2 solution is closer to 
those at the solution of the <(> 3 The result of solving 4>i (1 *" A> < 
is values that are on both 

<*><*-!>, and <t> k % 



Instead, the present invention defines all auxiliary func- 
tion 4^^ dimensional penalized cost function problem, lower 
order z* values and u* values determined previously by the 
reduction process. The *¥ k & k and its gradient is used for hill 
climbing to its peak. The z* values at the peak of the *¥ k into 
Problem Formulation [3.2] to determine the proper assign- 
ments. To define the explicitly makes the function ® k U* 
15 order sets of Lagrangian Coefficients with the expanded 
notation: & k U* U* +1 U* 4* is defined recursively, using the 



!3.6}(ai 

20 *3(" 3 ) = vz + X 4 " *^ ( " 3; M * uK) 



25 where v 2 is the solution value for the most recent 
2-dimensional penalized cost function problem. For k>3 



*h(« 3 



[3.61(b) 



30 



f <!>*(«*; u t+i If known 



35 



From the definition of & k k compared with Problem Formu- 
lation [3.21): 



40 



t k =o 



45 it follows that: 



*3(« 3 ')= V2+ 24 ^(W 3 ;"*. - : . ^(3?) 
'3=0 
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and [3.4] were solved at all possible values of u* the function an d w | tn mat Equation 3.5 is extended to: 

than the v t surface except at the peak as described in 
Equation 3.5 and its maximum occurs where the v* surface 
is minimum. Because the <Mc the <$> k & k v k surface. The & k 
minimization problem described by 13.11 (a). Let z* be the 
unknown solution to Problem Formulation 13.1] and note 
that: 



H'jtCu 3 , . . . , u*- 1 ;u*)£* s (tt*;tt* rt u^v^^u*) [3.7] 



C>*(u k )Svj k (z, 4 )Sv t _ l <z*) 



13.5] 



Consequently, the z* values at the peak of <J>* greatest lower 
bound on the cost of the relaxed problem), can be substituted 
into the k-dimensional penalized cost function to detennine 
the proper assignments. Consequently, the present invention 
attempts to find the maximum on the <$> k surface. However, 



This relationship means that either & k *¥ k hill climbing to 
60 update the Lagrangian Coefficients u. the preferred 
choice, however it is only available when the solution to 
Problem Formulation |3.2] is a feasible solution to Problem 
Formulation [3.1 1 (as in hill climbing from the second to 
third dimension which is why prior art worked). For sim- 
65 plicity in jimplementation, the function *¥ k that it equals 0* 
therefore always used, even for hill climbing from the 
second dimension. The use of the ¥ 
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The modification to Wolfe's algorithm uses the information 
generated for % ¥ k U* 3 U^~ 2 U"" 1 for calculating the sub- 
gradients. The definition of a subgradient v of U* 3 U** 
subdifferential set defined as: 

5 sy^iwR^noFidJ* 8 , . - - , y^Huysvu* 3 

U^U*))£v T (lT-U*) VUdT n> 

(where v T is the transpose of v) Next, a subgradient vector 
is determined from mis function. If t is the solution of 
10 Problem Formulation |3.2], then differentiating *¥ k U 3 U*" 1 
U* U k lk evaluating the result with respect to the current 
selection matrix z k yields a subgradient vector: 



15 



8 = 



climbing/peak detection are to determine the gradient of the 
x ¥ k and then move up the ¥ £ increasing portion of the 
gradient As shown in FIG. 6. any upward step on the *¥ k the 
4>„ U* a a is closer to the ideal set of u* values which 

correspond to the minimum of the v k function. While it is The iterative nature of the solution process at each dimen- 
possible to iteratively step up this x i f k and then determine if 20 sion yields a set of such subgradients. Except for a situation 
the peak has been reached, the present invention optimizes described below where the resultant averaged gradient does 
this process by determining the single step size from the not point toward the peak, the most recent set of such 
ing point at the minimum of ^ that will jump to the subgradients are saved and used as the "bundle" for the peak 
^o^anUthen calculating the u* values at the peak. Once the finding process for this dimension. For example, at the k 
^"u" s uM the peak at ¥ fc determined then the u* values can 25 level there is a bundle of subgradients of the *¥ k the mini- 
vfae^stittited into Problem Formulation [3.2] to determine mum of the ^ surface determined as a result of solving 
the proper assignment. (However, in a more realistic Problem Formulation [3.1] at all lower levels. This bundle 
example, where k is much greater than three, then the u* can be averaged to approximate the gradient. Alternately, the 
values at the peak of the ¥ along with ihe lower order u* previous bundle can be discarded so as to use the new value 
values and those assigned and yielded by the reduction steps 30 to initiate a new bundle. This choice provides a way to adjust 
are used to define the next higher level the process to differing classes of problems, i.e., when data 

is being derived from two sensors and the relaxation pro- 
¥ ceeds from data derived from one sensor to the other then the 

prior relaxation data for the first sensor could be detrimental 
35 to performance on the second sensors data. 
L2.3. S tep Size and Termination Criterion 
After the average gradient of the ¥ k determined, the next 
% ¥ t u k step is to determine a single step that will jump to the peak 

of the *¥ k is to first specify an arbitrary step size, and then 
40 calculate the value of V k gradient. If the *¥ k this probably 
means that the step has not yet reached the peak, 
determine the gradient of each ¥ Consequently, the step size is doubled and a new >¥ k value 

determined and compared to the previous one. This process 
is repeated until the new x ¥ k previous one. At that time, the 
45 step has gone too far and crossed over the peak. 
¥ Consequently, the last doubhng is rolled-back and that step 

size is used for the iteration. If the initial estimated *¥ k step 
* size is decreased by 50%. If this still results in a smaller ¥ fc 

previous step size is decreased by 25%. The following is a 
50 more detailed description of this process. 

With a suitable bundle of subgradients determined as just 
described. Wolfe's algorithm can be used to determine the 
effective subgradient (d) and the Upgraded value U^. 
From the previous iteration, or from an initial condition, 
"bundle") are determined at several points on the ¥ fc in the 55 ^ m exists a step length value (t)> ^ value< 
region of the starting point (Le.. minimum of § k ) and then 
averaged. Statistically* the averaged result should point u^/rtd 

toward the peak of the 4** Subgradient Algorithm (Wolfe75, t 
Wolfe79) for rninimization was previously known in another is calculated as an estimate of u ^ v To determine if the 
environment to determine a gradient of a non-smooth sur- 60 current step size is valid the we evaluate *¥ k u u - u + If the 
face using multiple subgradients and can be used with result represents an improvement then we double the step 
modification in the present invention. Wolfe's algorithm is size. Otherwise we halve the step size. In either case a new 
further described in **A Method of Conjugate Subgradients u + is calculated. The doubling or halving continues until the 
for MimnizingNondifferentiable Functions" page 147-173 step becomes too large to improve the result or until it 
published by Mathematical Programming Study 3 in 1975 65 becomes small enough to not degrade the result The result- 
and "Finding the Nearest Point in a Polytope" page 128-149 ing suitable step size is saved with d as part of the subgra- 
published by Mathematical Programming Study 11 in 1976. dient bundle. The last acceptable u + is assigned to u ^ 
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for / = 0,..., 
and 7 = 0, .... Ml 



|3.8| 



If the resulting 2-dimensional assignment problem. 
Minimize: ^2^^ 

Subject to: 7=1 M> 

/=o 

** 

7-1 M) 

^e{0, 1) 7=0 Afc,/ = 0,.„, 



10 



Three distinct criterion are used to determine when u*, is 
close enough to u k : 

1. The Wolfe's algorithm criterion of d=0 given that the 
test has been repeated with the bundle containing only 
the most recent subgradient. 5 

2. The difference between the lower bound <b k u* upper 
bound v k (z*, u A k) being less than a preset relative 
threshold. (Six percent was found to be an effective 
threshold for radar tracking problems.) 

3. An iteration count being exceeded. 
The use of limits on the iteration are particularly effective 

for iterations at the level 3 through n-1 levels, as these 
iterations will be repeated so as to resolve higher order 
coefficient sets. With limited iterations the process is in 
general robust enough to improve the estimate of upper 15 
order Lagrangian Coefficients. By limiting the iteration 
counts then the total processing time for the algorithm 
becomes deterministic. That characteristic means the pro- 
cess can be effective for real time problems such as radar, 
where the results of the last scan of the environment must be 20 
processed prior to the next scan being received. 
L2.4. Solution Recovery 

The following process determines if the u* values at what 
is believed to be the peak of the approximate the u* 
values at the peak of the corresponding <t> k function. This is 25 
done by determining if the corresponding penalized cost 
function is irrinimized. Thus, the u* values at the peak of the 
*¥ k Formulation |3.2] to determine a set of z assignments for 
k-1 points. During the foregoing reduction process. Problem 
Formulation |3.4| yielded a tentative list of k selected z 30 
points that can be described by their indices as; {(F,. . . 

where N 0 is the number of cost elements 
selected into the solution. One possible solution of Problem 
Formulation [3.1] is the solution of Problem Formulation 
[3.2] which is described as {(¥ t . . . ? w m* v . . 35 
with m* ( . . . ijt_ A as it was defined in Problem Formulation 
[3.3]. If this solution satisfies the k-th constraint set then it 
is the optimal solution. 

However, if the solution for Problem Formulation [3.2] is 
not feasible (decision 355), then the following adjustment 40 
process determines if a solution exists which satisfies the 
k-th constraint while retaining the assignments made in 
solving Problem Formulation [3.4J. To do this, a 
2-dimensional cost matrix is defined based upon all obser- 
vations from the k-th set which could be used to extend the 45 
relaxed solution. 
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has a feasible solution then the indices of that solution map 
to the solution of Problem Formulation [3.1] for the 



65 



k-dimensional problem. The first index in each resultant 
solution entry is the pointer back to an element of the {(F,. 
. . F^m* ( . . . ijt_i)}%:i list. That element supplies the first 
k-1 indices of the solution. The second index of the solution 
to the recovery problem is the k"' index of the solution. 
Together these indexes specify the values of z k that solve 
Problem Formulation |3.1| at the )s? h level. 

If Problem Formulation |3.9| does not have a feasible 
solution then the value of uki which was thought to represent 
u k p is noi representative of the actual peak and further 
iteration at the k"' level is required. This decision represent 
the branch path from step 214 and equivalent steps. 

1.3. Partitioning 

A partitioning process is used to divide the cost matrix 
that results from the Scoring step 154 into as many inde- 
pendent problems as possible prior to beginning the relax- 
ation solution. The partitioning process is included with the 
Problem Formulation Step 310. The result of partitioning is 
a set of problems to be solved, i.e., there will be p A problems 
that consist of a single hypothesis, p 2 problems that consist 
of two hypothesis, etc. Each such problem is a group in that 
one or more observations or tracks are shared between 
members of the group. 

In partitioning to groups no consideration is given to the 
actual cost values. The analysis depends strictly on the basis 
of two or more cost elements sharing the same specific axis 
of the cost matrix. In a 2-dimensional case two cost elements 
must be in the same group if they share a row or a column. 
If the two (elements are in the same row, then each other cost 
clement that is also in the same row, as well as any cost 
elements that are in columns occupied by members of the 
row must be included in the group. The process continues 
recursively. In literature it is referred to as "Constructing a 
Spanning Forest." The k-dimensional case is analogous to 
the 2-dimensional case. The specific method we have incor- 
porated is a depth first search, presented by Aho, Hopecraft, 
and Uilman, in "Design and Analysis of Computer 
Algorithms " section 5.2, published by Addison-Westley, 
Mass. 

The result of partitioning at level M is the set of problems 

described as ^P (> li=l p y and j=l, N}, where N is 

the total number of hypothesis. The problems labeled {P 4l li= 

1 pi \ are the cases where their is only one choice for 

the next observation at each scan and that observation could 
be used for no other track, i.e., it is a single isolated track. 
The hypothesis must be in included in the solution set and 
no further processing is required. 

As hypothesis are constructed the first clement is used to 
refer to the track id. Any problem in the partitioned set which 
does not have shared costs in the first scan represent a case 
where a track could be extended in several ways but none of 
the ways share an observation with any other track. The 
required solution hypothesis for this case is the particular 
hypothesis with the maximum likelihood. For this case all 
likelihoods are determined as was described in Scoring and 
the maximum is selected. 

In addition to partitioning at the M level, partitioning is 
applied to each subsequent level M-l, . . . , 2. For each 
problem that was not elirninated by either by the prior 
selection partitioning is repeated, ignoring observations that 
are shared in the last set of observations. Partitioning rec- 
ognizes that relaxation will eliminate the last constraint set 
and thus partitioning is feasible for the lower level problems 
that will result from relaxation. This process is repeated for 
all levels down to k=3. The full set of partitionings can be 
performed in the Problem Formulation Step 310, prior to 
initiating the actual relaxation steps. The actual 
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2-dimensional solver used in step 206 includes an equivalent 
process so no advantage would be gained by partitioning at 
the k=2 level. 

There are two possible solution methods for the re mainin g 
problems. "Branch and Bound" as was previously described, 5 
or the relaxation method that this invention describes. If any 
of the partitions have 5-10 possible tracks and less than 50 
to 20 hypotheses, then the prior art "Branch and Bound" 
algorithm generally executes faster than does the relaxation 
due to its reduced level of startup overhead. The "Branch 10 
and Bound" algorithm is executed against all re mainin g M 
level problems that satisfy the size constraint For the 
remainder the Relaxation algorithm is used. The scheduling 
done in Problem Formulation allows each Problem Formu- 
lation [3.2| cost matrix resulting from the first step of is 
relaxation to be partitioned. The resulting partitions can be 
solved by any of the four methods* isolated track direct 
inclusion, isolated track tree evaluation, small group 
"Branch and Bound" or an additional stage of relaxation as 
has been fully described. 20 

The partitioning after each Level of Lagrangian Relax- 
ation is effective because when a problem is relaxed, the 
constraint that requires each point to be assigned to only one 
track is eliminated (for one image at a time) . Therefore, two 
tracks previously linked by contending for one point will be 25 
unlinked by the relaxation which permits the point to be 
assigned to both tracks. The fruitfulness of partitioning 
increases for lower and lower levels of relaxation. 

The following is a more detailed description of the 
partitioning method Its application at all but the M level 30 
depends upon the relaxation process described in this inven- 
tion. The recursive partitioning is therefore a distinct part of 
this invention. The advantage of this method is greatly 
enhanced by the sparse cost matrix resulting from tracking 
problems. However the sparse nature of the problem 35 
requires special storage and search techniques. 

A hypothesis is the combination of the cost element c„ the 
selection variable z n and all observations that made up the 
potential track extension. It can be written as. h„={c„. z rt . 

{o n =okiik=l Mie{l Nk}}}. i.e- cost selection 40 

variable and observations in the hypothesis. ne{ 1 N} 

where N is the total number of hypotheses in the problem. 
While the cost and assignment matrices were previously 
referenced, these matrices are not effective storage mecha- 
nisms for tracking applications. Instead the list of all hypoth- 45 
esis and sets of lists for each dimension that reference the 
hypothesis set are stored. The hypothetical set in list form is: 

For each dimension k=l . . . M there exists a set of lists, with 50 
each list element being a pointer to a particular hypothesis: 

UHVkpri N^r^^.^ 

where is number of hypothesis containing the i-th 55 
observation from scan k. This structure is a multiply linked 
list in that any observation is associated with a set of pointer 
to all hypothesis it participates in, and any hypothesis has a 
set of pointers to all observations that formed it. (These 
pointers can be implemented as true pointers or indices 60 
depending upon the particular programming language 
utilized.) 

Given this description of the matrix storage technique 
then the partitioning technique is as follows: Mark the first 
list L kt and foDow out all pointers in that list to the indicated 65 

hypothesis for i=l Nki. Mark all located 

hypothesis, and for each follow pointers back the particular 
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L 0l for k=l M. Those Us if not previously marked get 

marked and. also followed out to hypothesis elements and 
back to Us. When all such Us or h's being located are 
marked, then an isolated sub-problem has been identified. 
All marked elements can be removed from the lists and 
stored as a unique problem to be solved. The partitioning 
problem then continues by again starting at the first residual 
L set. When none remain, the original problem is completely 
partitioned. 

Isolated problems can result from one source track having 
multiple possible extensions or from a set of source tracks 
contending for some observations. Because one of the 
indices of k (in our implementation it is k=l) indicate the 
source track then it is possible to categorize the two problem 
types by observing if the isolated problem includes more 
than one L-list from the index level associated with tracks. 
L4. Subsequent Track Prediction 
As noted above, after the points are assigned to the 
respective Hacks, some action is taken such as controlling 
take-offs and landings to avoid collision, advising airborne 
aircraft to change course, warning airborne aircraft of an 
adjacent aircraft in a commercial environment, or aiming 
and firing an anti-aircraft (or anti-missile) missile, rocket or 
projectile, or taking evasive action in a military environ- 
ment Also, the tracking can be used to position a robot to 
work on an object. For some or all of these applications, 
usually the tracks which have just been identified are 
extended or extrapolated to predict a subsequent position of 
the aircraft or other object. The extrapolation can be done in 
a variety of prior art ways; for example, straight line 
extensions of the most likely track extension hypothesis, 
parametric quadratic extension of the x versus time, y versus 
time, etc. functions that are the result of the filtering process 
described earlier, least square path fitting or Kalman filtering 
of just the selected hypothesis. 

The process of gating, filtering, and gate generation as 
they were previously described require that a curve be fit 
through observations such that fit of observations to the 
likely path can be scored and further so that the hypothetical 
target future location can be calculated for the next stage of 
gating. In tht implementation existing, quadratic functions 
have been fit through the measurements that occur within the 
window. A set of quadratics is used, one per sensor mea- 
surement To calculate intercept locations these quadratics 
can be converted directly into path functions. Intercept times 
are then calculated by prior methods based upon simulta- 
neous solution of path equations. 

The use of the fitted quadratics is not as precise as more 
conventional filters like the Kalman Filter. They are however 
much faster and sufficiently accurate for the relative scores 
required for the Assignment problem. When better location 
prediction is required then the assignment problem is 
executed to select the solution hypothesis and based upon 
the observations in those hypothesis the more extensive 
Kalman filter is executed. The result is tremendous compu- 
tation savings when compared with the Kalman Filter being 
run on all hypothetical tracks. 

Based on the foregoing, apparatus and methods have been 
disclosed for tracking objects. However, numerous modifi- 
cations and substitutions can be made without deviating 
from the scope of the present invention. For example, the 
foregoing ¥ u* values are used to eliminate the higher than 
approximation characteristic of the <f> k 
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be partitioned into tracks and false alarms. In the formula- 
¥ tion of track extensions a moving window over time of 

observations sets is used. The observation sets will be 
measurements which are: (a) assigned to existing tracks, (b) 
4> k 5 designated as false measurements, or (c) used for initiating 

tracks. However, note that alternative data objects instead of 
observation sets may also be fused such as in sensor level 
be as efficient it seems likely that the method would con- tracking wherein each senS or forms tracks from its own 
verge. Therefore the invention has been disclosed by way of measurements and tnen ^ Uacks from the sensors m fused 
example and not limitation, and reference should be made to io ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ 

the following claims to determine the scope of the present d both embodiments of the present invention ^ 

invention. , 1 1 , r . ^ * a + * • 

be used for this type of data fusion. 

E. AN ALTERNATIVE EMBODIMENT OF THE The data assignment problem considered presently is 

MULTI-DIMENSIONAL ASSIGNMENT 15 represented as a case of set partitioning defined in the 

SOLVING PROCESS following way. First, for notational convenience in repre- 

The following description provides an alternative second senting tracks, a zero index is added to each of the index sets 
embodiment of the multi-dimensional assignment solving in 1 5.2 1, a gap filler Z k 0 is added to each of the data sets Z(k) 

process of section LI. However, in order to clearly describe in I5.2|, and a "track of data" is defined as (z (i z tf ) 

the present alternative embodiment, further discussion is 20 where i k and z* tjt can now assume the values of 0 and z* 0 , 
first presented regarding the data assignment problem of respectively. A partition of the data will refer to a collection 
partitioning measurements into tracks and false alarms and of tracks of data or track extensions wherein each observa- 
in particular, the representation of multi-dimensional assign- tion occurs exactly once in one of the tracks of data and such 
ment problems. that all data is used up. Note that the occurrence of the gap 

ILL Formulation of the Assignment Problem The goal of 25 filler is unrestricted. The gap filler z* G serves several pur- 
this section is to explain the formulation of the data asso- poses in the representation of missing data, false 
ciation problems, and more particularly multi-dimensional observations, initiating tracks, and terminating tracks. Note 
assignment problems, that govern large classes of problems that a "reference partition" may be defined which is a 
in centralized or hybrid centralized-sensor level multisensor/ 3Q partition in which all observations are declared to be false, 
multitarget tracking. The presentation is brief; technical Next undef appropriate indep endence assumptions it can 

details are presented for both track initiation and mainte- . „. tUnt . 

. a w w tl . ... ^ £ be shown that: 

nance in A, B. Poore, Multi-dimensional assignment formu- 
lation of data association problems arising from multitarget 

tracking and multisensor data fusion, Computational Opti- P(T = y\Z ) ^ ^ ^ Li 

mizationandApplications,3(1994),pp.27-57,fornonma- f\r=Y°\z M ) oi ^r 

neuvering targets and in Ingemar J. Cox, Pierre Hansen, and 
Beta Julesz, eds., Partitioning Data Sets, DIMACS Series in 

Discrete Mathematics and Theoretical Computer Science. where L* (i i M is a likelihood ratio containing probabilities 
American Mathematical Society, Providence, R.L, v. Feb. ^ for detection, maneuvers, and termination as well as prob- 
19, 1995, pp. 169-198 for maneuvering targets. Note that the ability density functions for measurement errors, track ini- 
present formulation can also be modified to include target tiation and termination. Then if c*, i J ^-\nL k ll fA/ , 
features (e.g., size and type) into the scoring process 154. 1 * " 

The data assignment problems for multisensor and mul- \i\ Y \z M )\ v I 54 l 

titarget tracking considered here are generally posed as that 45 "H^z" ^ ^ Cii " iN ' 

of maximizing the posterior probability of the surveillance t'l.-^^ 
region (given the data) according to 



Maximize 



\ 1511 r x 



Using [5.3 ] and the zero-one variable z*^ ^1 if (i^ . . 



50 



e y 



|5.5| 



where Z M represents M data sets, y of the data (and thus 
induces a partition of the data), T* the finite collection of all 
such partitions, V 55 
* 

Y(ryZ M ) is the posterior probability of a partition y true y iP £ „ £ Cjj Jtf?ii JM 

given the data Z M . The term partition is defined below. i t =o * w «o 

Consider M observation sets Z(k) (k=L . * . , N) each of 
N t observations {z*^}"*^, and let Z M denote the cumulative eo y V ... V ^ ^=i t l ^ 

data set defined by * 2 =o i M «o 1 

(o y ... x v ... x zi t 

respectively. In multisensor data fusion and multitarget , 1= g 1=0^+1=0 ^=0 

tracking the data sets Z(k) may represent different classes of 65 
objects, and each data set can arise from different sensors. 
For track initiation the objects are measurements that must 



for i t = l N k and k = 2, M - 1, 
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-continued 

'1=° 

€{0 > lf ' fora11 il iM 



where Co 0 is arbitrarily defined to be zero. Here, each 
group of sums in the constraints [5.5 J (b)-[5.5] (e) repre- 
sents the fact that each non-gap filler observation occurs 
exactly once in a i4 track of data." One can modify this 
formulation to include multiassignments of one, some, or all 
the actual observations. The assignment problem [5.5] is 
changed accordingly. For example, if z k lk is to be assigned no 
more than, exactly, or no less than n* u times, then the "=1" 
in the appropriate one of the constraints |5.5[ (b)-[5.5] (d) 
is changed to =, =n* respectively. 

Expressions for the likelihood ratios L* - - « M such as \ 5.5] 
are well known. In particular, discussions of these ratios may 
be found in: A.B. Poore, Multi-dimensional assignment 
formulation of data association problems arising from mul- 
titarget tracking and multisensor data fusion. Computational 
Optimization and Applications, 3 (1994), pp. 27-57; and 
Ingemar J. Cox, Pierre Hansen, and Beta Julesz, eds.. 
Partitioning Data Sets, DIMACS Series in Discrete Math- 
ematics and Theoretical Computer Science. American Math- 
ematical Society, Providence, R.L. v. 19, 1995. pp. 169-198. 
Furthermore, the likelihood ratios are easily modified to 
include target features and to account for different sensor 
types. Also note that in track initiation, the M observation 
sets provide observations from M sensors, possibly all the 
same. Additionally note that for track maintenance, a sliding 
window of M observation sets and one data set containing 
established tracks may be used. In this latter case, the 
formulation is the same as above except that the dimension 
of the assignment problem is now M+l. 

n.2. Overview of the New Lagrangian Relaxation Algo- 
rithms 

Having formulated an M-dimensional assignment prob- 
lem [5.5], we now turn to a description of the Lagrangian 
relaxation algorithm. Subsection IL2.1 below presents many 
of the relaxation properties associated with the relaxation of 
an ifdimensional assignment problem to an m-dimensional 
one via a Lagrangian relaxation of n-m sets of constraints 
wherein m<n<M and preferably in the present invention 
embodiment c-m>l. Although any n-m constraint sets can 
be relaxed, the description here is based on relaxing the last 
n-m sets of constraints and keeping the first m sets. Given 
either an optimal or suboptimal solution of the relaxed 
problem, a technique for recovering a feasible solution of the 
n-dimensional problem is presented in subsection IL2.2 
below and an overview of the Lagrangian relaxation algo- 
rithm is given in subsection IL2.3. 

The following notation will be used throughout the 
remainder of the work. Let M be an integer such that M^3 
and let ne{3, , . . M}. The n-dimensional assignment 
problem is 
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-continued 

for i t = 1 N k 

and A = 2, ...»« - 1, 



, € {0, 11 



for all i x . 



15 To ensure that a feasible solution of [6.1] always exists, all 
variables with exactly one nonzero index (i.e.. variables of 
the form iz^o 0tt0 0 for i k #0) are assumed free to be 
assigned and the corresponding cost coefficients are well- 
defined. 

20 

H 2.1. The Lagrangian Relaxed Assignment Problem 

The n-dimensional assignment problem [6.1] has n sets of 
constraints. A (N^+l)-dimensional multiplier vector associ- 

25 ated with the k-th constraint set will be denoted by u*=(u* 0 , 

u*j u k M f with u* o =0 and k=l n. The 

n-dimensional assignment problem [6.1] is relaxed to an 
m-dimensional assignment problem by incorporating n-m of 
the n sets of constraints into the objective function [6.1] (a). 

30 Although any constraint set can be relaxed, the description 
of the relaxation procedure for [6.1] will be based on the 
relaxation of the last n-m sets of constraints. The relaxed 
problem is 
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1 = Minimize 



[6.2] 



N k-l N k+l 



k=m+l i t =0 



E-2h-'\M< 

A N m 



Subject to: 



1 N u 



i!-0 ^=0,^-0 

for k = 1 and k = 2, ™, 

€10.1} for aD h i„. 



(a) Minimize: v„(z) = X""Z < ?i™t.*i™< 1 . 

<b) Subjectto: £ "* = U h = 1 



60 wherein we have added the multiplier u* o =0 for notational 
* 6 convenience. Thus, the multiplier R with u* o =0 is fixed. 

An optimal (or suboptimal) solution of [6.2] can be 
constructed from that of an m-dimensional assignment prob- 
65 lem. To show this, def ine f or each (i x i^) an index 

Om^l J«) KW(*I un> Uh a 

new cost function c m i{ ^ m hy 
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-continued 

Um+i jw)=argmin I 6 - 3 ! h^.,.^ = 0 if ... i„) * — jb) and 



(jt =0, .... tf t and Jt = m+ 1, 



iii UMO,... f O) 

* 



4-. < .-«f,.- Wwl _i+ Z for(i ' w*«a-.«) 



10 



* =0 " Then w" is a feasible solution of the Lagrangian relaxed 

problem |6.21 and 

(If (j^j j w n) is not unique, choose the smallest such 15 

im+i* amongst those (n-m)-tuples with the same j m+1 choose ^ -f m + l n ^ w";tf" +1 , . . . , uv^i^^, 

the smallest j 2 etc so that Q i gJ*.™^* If. in addition. V" is optimal for |6.4], then W is an optimal 

defined.) Using the cost coefficients denned in this way, the so i ut ion Q f 15 21 and 
following m-dimensional assignment problem is obtained: 



With the exception of one equality being converted to an 
A inequality, the following Fact is a converse of Fact A.2 and 

S Aj " E ^ *'* tm ^ "' im nas ^ so ^ cen shown to be true. 

'i=o irt,=0 25 Fact A3. Let w" be a feasible solution to problem 16.2] 

and define W m by 



Subject to: 



"1 



30 



^ — ^ /„=(> 



Z-ZZ •£>«--'• 

ii-o i*-i-o** + i-o ^ o = 0 if {iu ...j m ) = ( 0 0) and 



4».<W..t + E <>Oforall(w,...,U 
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for i t = 1, JV t and Jfc = 2, JV„ 

> - E «s— 

^ '»-i=° <„.o = 1 if («i W = (0, .... 0) and 

zT...( m €{0, 1} for all i x t m . A 4 

i=m+I 

40 

As an aside, observe that any feasible solution z" of [6.1] 

yields a feasible solution z m of [6.4] via the construction Then w 7 " is a feasible solution of the problem [6.4] and 

{1 if 2?,...^ 4, = 1 for some (i M+l i n ) 
45 If, in addition, w" is optimal for [6.2], then w™ is an optimal 
0 solution of [6.4], m+1 " 

Uw";^, u-M^.Z^^, and <J> W " 

Thus, the mn-dlmensional assignment problem [6.41 has at 

least as many feasible solutions of the constraints as the 5Q l > * • • > ^ >-^j™i2"V»V 

original problem |6.1 1. n 2 2 irhe Rccovery procedure 

The following Fact A.2 has been shown to be true. It states next objective is to explain a recovery procedure or 

matanoptimalsolutionofProblemFormulation|6.21canbe me thod for recovering a solution to the n-drmensional 
computed from that of Problem Formulation [6.41 problem of [6.1 1 from a relaxed problem having potentially 
Furthermore, the converse to this fact is provided in Fact 55 substantially fewer dimensions than [6.1J. Note that this 
A3. Moreover, if the solution of either of these two prob- aspect of me alternative embodiment of the present inven- 
iems [6.2 1 or |6.4| is e-optimal (i.e., the objective function tion is substantially different from the method disclosed in 
associated with the suboptimal solution is within "e", of the the ^ embodiment of the multi-dimensional assignment 
objective function associated with the optimal solution), solving process of section 1. 1 of this specification. Further, 
then so is the other. 60 this alternative embodiment provides substantial benefits in 

Fact A.2. Let \f be a feasible solution to problem [6.4] terms of computational efficiency and accuracy over the first 
and define w" by embodiment, as will be discussed. Thus, given a feasible 

(optimal or suboptimal) solution w™ of [6.4] (or w" if 
*? i = im if <*»+i. — 4> = Uw+i» — In) and l 6 5 3 Problem Formulation [6.2] is constructed via Fact A.2), an 

65 explanation is provided here regarding how to generate a 

{ix im) * (0 * *" * 0) feasible solution z" of [6. 1] which is close to w" 1 in a sense 

to be specified. We first assume that no variables in [6. 1] are 
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preassigned to zero (this assumption will be removed 
shortly). The difficulty with the solution w" is that it need not 
satisfy the last n-m sets of constraints in [6.1]. Note, 
however, that if w m is an optimal solution for [6.4] and w" 
(constructed as in Fact A.2) satisfies the relaxed constraints, 
men w n is optimal for |6.1]. The recovery procedure 
described here is designed to preserve the 0-1 character of 
the solution of [6.4] as far as possible. That is. if w m Ij . . . 

i m =l and ilw*0 for at least one 1=1 m. then the 

corresponding feasible solution z" of [6.1] is constructed so 
that z",.. . . i n =l for some (w- . . i n ). Note that by mis 



reasoning, variables of the form z" 0 



_ in can be 

assigned to a value of 1 in the recovery problem only if 
w» 0 0 =1. However, variables z" 0 0w . mlm will be treated 
differently in the recovery procedure in that they can be 
assigned 0 or 1 independent of the value v/" 0 , . <>. This 
increases the feasible set of the recovery problem, leading to 
a potentially better solution. 

Let {(?,. . . . iV \ / 
w (or the first m indices of wl constructed as in Fact A.2) 
such that w"V . . and (^i ^J*^ • • • > °) • Set 



,i°J=<0 0) for j=0 and define 



for ijt = 



= 0 N t \ 

1,..., n;j = Q, . 



16-71 



, No- 



Let Y denote the solution 
assignment problem: 



Minimize 



Subject to 



-m+l 



"0 "m+l » 



[6.8] 



= 1» 



for ijt = 
and k-- 



m+ L 



*1- 



1 if d'i ft) for some 

;-0,...,Afe and ^ r ..^ = l 



0 otherwise. 
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zero. To ensure that a feasible solution exists, we now only 



need ensure that the variables z" 



are free for. 



10 



15 



j=(U N 0 . (Recall that all variables of the form 

z rt 0 l o are free (to be assigned) and all coefficients of the 

form c^o " . ... 0 arc well defined for k=l n.) This 

is accomplished as follows. If the cost coefficient c"^ . . . 
¥„ 0 ... 0 is well defined and z"JJ / m 0 ... 0 is free, then 
define , . epe"^ • • / J) . . ; 6 with z^ 1 * . 0 being 

free. Otherwise, since all variables of the form z 0 ijt . 0 
are known feasible and have well-defined costs, put 



De an enumeration of indices of 



20 



25 



of the (n-m+l)-dimensional 



30 



35 



40 



45 



^ {0, 1} for aU h- 

The recovered feasible solution z" of [6.1] corresponding to 
the multiplier set {u m u rt } is then defined by 



50 



[6.9] 



55 



60 



n-m+1 



= £ 



This recovery procedure is valid as long as all cost coeffi- 
cients c" are defined and all zero-one variables in z" are free 
to be assigned. Note that modifications are necessary for 
sparse problems. If the cost coefficient c\ . J mi ^ v ^ is well 
defined and the zero-one variable z\ . 'J^*. . . ^ is free 

, then define d % - m * 1 ii ... . . /„= 65 



to be assigned to zero or one. 



ifmWl* 



. i„ as in [6.7] with z* 



to be assigned. Otherwise. z J 



-^ w ..i n being free 
t n is preassigned to 



11.3. The Multi-Dimensional Lagrangian Relaxation 
Algorithm for the n-Dimensional Assignment Problem 

Starting with the M-dimensional assignment problem 
16.1], i-e.. n=M. the algorithm described below is recursive 
in that the M-dimensional assignment problem is relaxed to 
an m-dimensional problem by incorporating (n-m) sets of 
constrains into the objective function using Lagrangian 
relaxation of this set. This problem is maximized with 
respect to the Lagrange multipliers, and a good suboptimal 
solution to the original problem is recovered using an 
(n-m+1) -dimensional assignment problem. Each of these 
two (the m-dimensional and the (n-m+1) -dimensional 
assignment problems) can be solved in a similar manner. 

More precisely, reference is made to FIG. 8 which pre- 
sents a flowchart of a procedure embodying the multi- 
dimensional relaxation algorithm, referred to immediately 
above. This procedure, denoted MULTE i3 DIM__RELAX in 
FIG. 8, has three formal parameters, n, PROB_ 
FORMULATION and ASSIGNMT_SOLUTION, which 
are used to transfer information between recursive instan- 
tiations of title procedure. In particular, the parameter, n. is an 
input parameter supplying the dimension for the multi- 
dimensional assignment problem (as in [6.1]) which is to be 
solved (i.e.. to obtain an optimal or near-optimal solution). 
The parameter, PROB_FORMULATION, is also an input 
parameter supplying the data structure(s) used to represent 
the n-dimensional assignment problem to. be solved. Note 
that PROB__FORMULAIION at least provides access to the 
cost matrix., [c"], and the observation sets whose observa- 
tions are to be assigned. Additionally, the parameter. 
ASSIGNMT_SOLUTIOR is used as an output parameter 
for returning a solution of a lower dimensioned assignment 
problem to an instantiation of MULTI_DIM_RELAX 
which is processing a higher dimensioned assignment prob- 
lem. 

A description of FIG. 8 follows. Assuming MULTL_ 
DIM_RELAX is initially invoked with M as the value for 
the parameter n and the PROB_FORM ULATION having a 
data structure^) representing an M-dimensional assignment 
problem as; in [5.51. in step 500 an integer m, 2^m<n is 

chosen such that the constraint sets or dimensions m+l 

n are to t>e relaxed so that an m-dimensional problem 
formulation as in [6.21 results. In step 504 an initial approxi- 
mation is provided for {u"*"^ u"o56 . Subsequently, in 

step 508 the above initial values for {u""" 1 ^. . . ,u n 0 } are 

used in an iterative process in determining (u 17 ^ 1 u rt ) 

which rnaximizes {^Y"* 1 " *eR Wvn and k=m+l. . . . ji} for 
a feasible solution w^ subject to the constraints of Problem 
Formulation [6.2]. Note mat by maxiinizing 0„ m Y m ' 1 * 
constraints that were relaxed are being forced to be satisfied 
and in so doing information is built into a solution of this 
function for solving the input assignment problem in 
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"PROB_FORMULATION." Further note that a non- Thus, in step 544. the n-m+1 solution is used to solve the 
smooth 10 optimization technique is used here and that a n-dimensional assignment problem as discussed regarding 
preferred method of determining a maximum in step 508 is equations 16.9 J. Finally, in steps 548 and 552 the solution to 
the bundle-trust method of Schramm and Zowe as described the n-dimensional assignment problem is returned to the 
in H. Schramm and J. Zowe, A version of the bundle idea for 5 calling program so that for example, it may be used in 
minimizing a non-smooth function: Conceptual idea, con- taking one or more actions such as (a) sending a warning to 
vergence analysis, numerical results, SIAM Journal on aircraft or sea facility; (b) controlling air traffic; (c) control- 
Optimization, 2 (1992), pp. 121-152. This method, along ling anti-aircraft or anti-missile equipment; (d) taking eva- 
with various other methods for determining the maximum in sive action; or (e) surveilling an object, 
step 508, are discussed below. io H.3.L Comments on the Various Procedures Provided by 

Also note that for m>2, a solution to the optimization FIG. 8 
problem of step 508 is NP-hard and therefore cannot be There are many procedures described by FIG. 8. One such 
solved optimally. That is, there is no known computationally procedure is the first embodiment of the multi-dimensional 
tractable method for guarantying an optimal solution. Thus, assignment solving.process of section LI. That is, by speci- 
there are two possibilities: either (a) allow m to be greater 15 fying m=n -1 in step 500, a single set of constraints is 
than 2 and use auxiliary functions similar to those disclosed relaxed in step 508. Thus, one set of constraints incorporated 
in the first embodiment of the k-dimensional assignment is inta the objective function via the Lagrangian problem 
solver 300 in section I to compute a near-optimal solution, formulation, resulting in an m=n-l dimensional problem, 
or (b) make m=2 wherein algorithms such as the forward/ The relaxed problem is subsequently maximized in step 512 
reverse auction algorithm of D. P. Bertsekas and D. A. 20 with respect to the corresponding Lagrange multipliers and 
Castafion in their paper, A forward/reverse auction algorithm then a feasible solution is reconstructed for the 
for asymrmnetric assignment problems. Computational n-dimensional problem using a two-dimensional assignment 
Optimization and Applications, 1 (1992), pp. 277-298, problem. The second procedure provided by FIG. 8 is a 
provides an optimal solution. novel approach which is not suggested by the first embodi- 
If option (a) immediately above is chosen, then the 25 ment of section LI. In fact, the second procedure is some- 
auxiliary functions are used as approximation functions for what of a mirror image of the first embodiment in that n-2 
obtaining at least a near-optimal solution to the optimization sets of constraints are simultaneously relaxed, yielding 
problem of step 508. Note that the auxiliary functions used immediately an m=2 dimensional problem in step 512. Thus, 
depend on the value of m. Thus, auxiliary functions used a feasible solution to the n-dimensional problem is then 
when m=3 will likely be different from those for m=4. But 30 recovered using a recursively obtained solution to an n-1 
in any case, the c^tirnization procedure is guided by using dimensional problem via step 540. In this case, the function 
the merit function, O^y 3 " 2-dimensional assignment prob- values and subgradients of cD^y 3 " optimally via a two 
lem for guiding the maximization process. dimensional assignment problem. The significant advantage 
Alternatively, if option (b) above is chosen, then two here is. that there is no need for the merit or auxiliary 
important advantages result: the optimization problem of 35 functions, x Y k embodiment of the multi-dimensional assign- 
step 508 can be always solved optimally and without using ment solving process of section LI above and also there is 
auxiliary approximation functions. Thus, better solutions to no need for the more general merit or auxiliary functions 
the original M-dimensional assignment problem are likely 4*^ subsection II.4.2 below. Further, note that all function 
since there is no guarantee that when the non-smooth values and subgradients used in the non-smooth niaximiza- 
optimization techniques are applied to the auxiliary func- 40 tion process are computed exactly (i.e., optimally) in this 
tions the techniques will yield an optimal solution to step second procedure. Moreover, problem decomposition is now 
508. Furthermore, it is important to note that without aux- carried out for the n-dimensional problem; however, decom- 
iliary functions, the processing in step 508 is both concep- position of the n-1 dimensional recovery problem (and all 
tually easier to follow and more efficient. _ lower recovery problems) is performed only after the prob- 

Subsequently, in step 512 of HG. 8, the solution (u"*" 1 , 45 lem is formulated. 

. . . ♦ u)" and w" is used in determining an optimal solution, Between these two procedures are a host of different 

w™ to Problem Formulation |6.4 1 as generated according to relaxation schemes based on relaxing n-m sets of constraints 

[6.3] and Fact A.3. to an m-dimensional problem (2<m<n), but these all have 

In step 516, the solution w" is used in defining the cost the same difficulties as the procedure for the first embodi- 

matrix a*""* 1 as in |6.7]. Subsequently, if n-m+l=2, then 50 ment of section L 1 in that the relaxed problem is an NP-hard 

the assignment problem [6.8] may be solved straightfor- problem. To resolve this difficulty, we use an auxiliary or 

wardly using known techniques, such as forward/reverse merit function ^^Ynotation^^i^m^-l, the recovery 

auction algorithms. Following this, in step 528, the solution procedure is still based on an NP-hard n-m+1 dimensional 

to the 2-dimensional assignment problem is assigned to the assignment problem. Note that the partitioning techniques 

variable ASSIGNMT_SOLUTION and in step 532 55 similar to those discussed in Section 1.3 may be used to 

ASSIGNMT_SOLUTION is returned to a dimension three identify the assignment problem with a layered graph and 

level recursion of MULTI_DIM_RELAX for solving a then to find me disjoint components of this graph. In general, 

3-dimensional assignment problem. all relaxed problems can be decomposed prior to any non- 

Aiternatively, if in step 520, n-m+l>2, then in step 536 smooth computations because their structure stays fixed 

the data structure(s) representing a problem formulation as 60 throughout the-algorithm of FIG. 8. However, all recovery 

in [6.8] is generated and assigned to the parameter, PROB„ problems cannot be decomposed until they are formulated, 

FORMULATION. Subsequently, in step 540 a recursive as their structure changes as the solutions to the relaxed 

copy of MULH_DIM_RELAX is invoked to solve the problems change. 

lower dimensional assignment problem represented by E.4. Details, and Refinements Relating to the Flowchart 

PROB_FORMULAI10N. Upon the completion of step 65 of FIG. 8 

540, the parameter, ASSIGNMT _JSOLUTION, contains the Further explanation is provided here on how various steps 

solution to the n-m+1 dimensional assignment problem* of FIG. 8 are solved. Note that the refinements presented 
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here can significantly increase the speed of the relaxation 
procedure of step 508. 
n.4.1 Maximization of the Non-smooth Function 
One of the key steps in the Lagrangian relaxation algo- 
rithm in section IL3 is the solution of the problem 

Maximize {S^?* 1 n *R M *1; Jt=m+1, . . . , n } [7.1] 

where u o *=0 for all k=m+l. . . . jl The evaluation of 
Qrf** 1 "corresponding niirumization problem [6.2] Zn var- 
ies for each instance of (u"*" 1 u n ). The following 

discussion provides some properties of these functions. 

Fact A.4. Let u" 1 * 1 u" be multiplier vectors associated 

with the (m+l) w through the n** set of constraints [6.1], let 
®mn Y 11 " objective function value of the n-dimensional 
assignment problem in equation [6.1], let z n be any feasible 
solution of [6.1], and let z* be an optimal solution of [6.1]. 
Then, O^y m+i n {xf* 1 u n }and 



Furthermore. 



[7.2a] 



[7.2a] 



for all u^R^ with u o *=0 and k=m. . . . n. 

Most of the procedures for non-smooth optimization are 
based on generalized gradients called subgradients, given by 

the following definition. Definition. At uK"""* 1 u n ) the 

set SO^y subdifferential of 

Ir. . . xR M "U} [7.3] 

A vector qeSO^y 

If z" is an optimal solution of [6.2] computed during 
evaluation of ^^y £ yields the following i n -th com- 
ponent of a subgradient g" of for S^y 

«S-o [7.4] 
for ijt = 1, Alt ami k - m + l t n. 



If z" is the unique optimal solution of [6.2], SO^y 
unique, then there are finitely many such solutions, say 
z n (l), . . . ,z n (K). Given the corresponding subgradients, 

g 1 g* the subdifferential 6<1>y {g 1 g*}. 

II.4.2 Mathematical Description of a Merit or Auxiliary 
Function 

For real-time needs, one must address the fact that the 
non-smooth optimization problem of step 508 (FIG. 8) 
requires the solution of an NP-hard problem for m>2. One 
approach to this problem is to use the following merit or 
auxiliary function to decide whether a function value has 
increased or decreased sufficiently in the line search or trust 
region methods: 

S";** 1 «")= P*5] 

Qnmiu"* 1 **) if m = 2 

or (6.2) is solved optimally, 
^(fi 3 S*;jr" +1 ,...,0 otherwise. 



where the multipliers u 3 u m that appear in lower order 

relaxations used to construct (suboptimal) solutions of the 



40 



m-dimensiional relaxed problem (6.2) have been explicitly 
included. Note that 4^ y be solved optimally if m=2. For 
sufficiently small problems (6.2) or (6.4), one can more 
efficiently solve the NP-hard problem by branch and bound. 
5 This is the reason for the inclusion of the first case; 
otherwise, the relaxed function ® 2n me function 
provides a lower bound for the optimal solution follows 
directly from Fact A.4 and the following fact 



10 



Fact A.5. Given the definition of y 

vp 3 m m+1 yjM-1 n jp ^1 



1!" [7.6] 

The actual function value used in the optimization phase is 
%n« Y with the solution z r ". . . i n being a suboptimal solution 
constructed from a relaxation procedure applied to the 

15 m-dimensional problem. Again, the use of these lower order 
relaxed problems is the reason for the inclusion of the 
multipliers u 3 . . . , u m . 

To explain how the merit function is used, suppose there 
is a current multiplier set (u^^ 1 u ol /) and it is 

20 desirable to update to a new multiplier set (u„ m+1 



via (u 

n 



/Hu^r 1 cwa^ 



A n oid m oi?* 1 o«" A y a J old m & obtained during the 

relaxation process used to compute a high quality solution to 
the relaxed m-dimensional assignment problem (6.2) at 

25 (u* 1 * 1 u n )=(u £>w m+i u oi /). The decision to accept 

,) is then based on ^ 3 old m old u oi /*K 
"Astopping criteria 



(u/^ew. . . ,u" 



», 3 fTt fftrhl 

rani new new new new 



commonly used in line searches. Again, (u^ . . . . M MW m ) 
represents the multiplier set used in the lower level relax- 
30 ation procedure to construct a high quality feasible solution 
to the m-dimensional relaxed problem (6.2) at (u"""* 

""M^^ 1 U ne»")' Trie P 0iDt is that timt 0JQe 

changes (u: m_t " 1 u n ) and uses the merit function *F m y 3 

generally change the lower level multipliers (u 3 , 

35 u m ). 

An illustration of this merit function for m=n-l is given 
in "PARTITIONING MULTIPLE DATA SETS. MULTIDI- 
MENSIONAL ASSIGNMENTS AND LAGRANGIAN 
RELAXATION", by Aubrey B. Poore and Nenad Rijavec, 

40 and in "QUADRATIC ASSIGNMENT AND RELATED 
PROBLEMS. DIMACS SERIES IN DISCRETE MATH- 
EMATICS AND THEORETICAL COMPUTER 
SCIENCE ", in American Mathematical Society, Providence, 
R.L, Vol. 16, 1994, pp. 25-37 edited by Panos M. Pardalos 

45 and Henry Wolkowicz. 

n.4.3 Non-smooth (^timization Methods 
By Fact A.4 the function A affine. and concave, so 
that the negative of A Thus the problem of maximizing 
O^y A optimization. Since there is a large literature on such 

50 problems, only a brief discussion of the primary classes of 
methods for solving such problems is provided. A first class 
of methods, known as subgradient methods, are reviewed 
and analyzed in the book by N. Z. Shor, Minimization 
Methods for Non-Differentiable Functions* Springer- Veriag, 

55 New York. 1985. Despite their relative simplicity, subgra- 
dient methods have some drawbacks that make them inap- 
propriate for the tracking problem. They do not guarantee a 
descent direction at each iteration, they lack a clear stopping 
criterion, sind exhibit poor convergence (less than linear). 

60 A more recent and sophisticated class of methods are the 
bundle methods; excellent developments are presented in the 
books of Hiriart-Urruty and Lemarechal: J.-B. Hiriart- 
Urruty and C. Lemarechal, Convex Analysis and Minimiza- 
tion Algorithms I and II. Springer- Verlag, Berlin, 1993; and 

65 the book of K. C Kiwiel: Methods of Descent for Non- 
Differentiable Optimization, Lecture Notes in Mathematics^ 

Vol 1133, Springer-Verlag, Berlin, 1985. Bundle methods 
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retain a set (or bundle) of previously computed subgradients new observation set) arrives from a sensor, one can obtain an 

to determine the best possible descent direction at the initial primal feasible solution by matching new reports to 

current iteration. Because of the non-smoothness of <b mn existing tracks via a two dimensional assignment problem, 

may not provide a "sufficient" decrease in <l> wt bundle This is known as the track-while -scan (TWS) approach, 

algorithms take a "null" step, wherein the bundle is enriched 5 Thus, an initial primal solution exists and then we use the 

by a subgradient close to u*. As a result, bundle methods are above relaxation procedure is to make improvements to this 

non-ascent methods which satisfy the relaxed descent con- TWS solution. Also for track maintenance, multipliers are 

dition O rofl y Jk+1 A (J^y^A k+l ^u^ These methods have been available from a previously solved and closely related 

shown to have good convergence properties. In particular, multi-dimensional assignment problems for all but the new 

bundle method variants have been proven to converge in a 10 observation set. 

finite number of steps for piecewise affine convex function- n.4.6. Local Search Methods 

als in K.. Kiwiel, An Aggregate Subgradient Method for Given a feasible solution of the multi-dimensional assign- 
Non-smooth Convex Minimization, Mathematical Program- ment problem, one can consider local search procedures to 
rning 27, (1983) pp. 320-341. and H. Schramm and J. Zowe, improve this result, as described in C. H. Papadimitriou and 
A version of the bundle idea for minimizing a non-smooth 15 K. Steiglitz, Combinatorial Optimization: Algorithms and 
function: Conceptual idea, convergence analysis, numerical Complexity* Prentice-Hall, Englewood Cliffs, N.J., 1982, 
results, SIAM Journal on Optimization, 2 (1992), pp. and J. Pearl, Heuristics: intelligent search strategies for 
121-152. computer problem solving, Addison- Wesley, 1984. One 
All of the above non-smooth optimization methods method is the idea of k-interchanges. Since for sparse 
require the computation A eStf^y A the computation 20 problems only those active arcs that participate in the 
of the relaxed cost coefficients |6.3|. In both the first and solution are stored, it is difficult to efficiently identify 
second procedures discussed in section H3. 1, maximization eligible variable exchange sets with some data structures for 
of 4>2, t y A In the most efficient implementations presently solving the assignment problem. However, a local search 
known of these two procedures, it was found that at least that may be more promising is to use the two dimensional 
95% of the computational effort of the entire procedure is 25 assignment solver in the following way. Given a feasible 
spent in the evaluation of the relaxed cost coefficients [ 6.3 1 solution to the multi-dimensional assignment problem, the 
as part of computing <J> 2rt y A that makes as efficient use of the indices that correspond to active arcs in the solution are 
subgradients and function values as is possible, even at the enumerated. Subsequently, one coordinate is freed to for- 
cost sophisticated manipulation of the subgradients. In mulate a two dimensional assignment problem with one 
evaluating three different bundle procedures: (a) the conju- 30 index corresponding to the enumeration and the other to the 
gate subgradient method of Wolfe used in section LI of the freed coordinate, and then solve a two dimensional assign- 
first embodiment of the present invention; (b) the aggregate ment problem to improve the freed index position, 
subgradient method of Kiwiel (described in K. C. Kiwiel, An IL4.7. Problem Decomposition 

Aggregate Subgradient Method for Non-smooth Convex The procedures described thus far are all based on relax- 

Minimizjation* Mathematical Programming 27, (1983) pp. 35 ation. Due to the sparsity of assignment problems, however, 

320-341); and (c) the bundle trust method of Schramm and frequently a decomposition of the problem into a collection 

Zowe (described in H. Schramm and J. Zowe, A version of of disjoint components can be done wherein each of the 

the bundle idea for rninimizing a non-smooth function: components can be solved independently. Due to the setup 

Conceptual idea, convergence analysis, numerical results, costs of Lagrangian relaxation, a branch and bound proce- 

SIAM Journal on Optimization, 2 (1992), pp. 121-152), it 40 dure is generally more efficient for small components, say 

was determined that for a fixed number of non-smooth ten to twenty feasible 0-1 variables (i.e., z ti t J. Otherwise, 

iterations, say, ten, the bundle-trust method provides good the relaxation procedures described above are used. Perhaps 

quality solutions with the fewest number of function and the easiest way to view the decomposition method is to view 

subgradient evaluations of all the methods, and is therefore the reports or measurements as a layered graph. A vertex is 

the preferred method. 45 associated with each observation point, and an edge is 

n.4.4 The Two Dimensional Assignment Problem allowed to connect two vertices only if the two observations 

The forward/reverse auction algorithm of Bertsekas and belong to at least one feasible track of observations. Given 

Castanon (as described in D. P. Bertsekas and D. A. this graph, the decomposition problem can then be posed as 

Castafion, A forward/reverse auction algorithm for asym- that of identifying the connected subcomponents of a graph 

metric assignment problems. Computational Optimization 50 which can be accomplished by constructing a spanning 

and Applications, 1 (1992), pp. 277-298) is used to solve the forest via a depth first search algorithm, as described in A. 

many relaxed two dimensional problems that occur in the V. Aho, J. E. Hopcroft, and J. D. Ullman, The Design and 

course of execution. Analysis of Computer Algorithms, Addison- Wesley Publish- 

IL4.5. Initial Multipliers and Hot Starts ing Company, Reading, Mass., 1974. Additionally, decom- 

The effective use of "hot starts" is fundamental for 55 position of a different type might be based on the identiri- 

real-time applications. A good initial set of multipliers can cation of hi- and tri-connected components of a graph and 

significantly reduce the number of non-smooth iterations enumerating on the connections. Here is a technical expla- 

(and hence the number of 0^ A quality recovered solution. nation. Let 

Further, there are several ways that multipliers can be 2^z^...i rt iziii2...inisnotpreassignedtozcio} 
reused. First, if a n-dimensional problem is relaxed to a 60 

m-dimensional problem, relaxation provides the multiplier denote the set of assignable variables. Define an undirected 

set {u m+1 , u"^ 2 , u n }. These can be used as the initial graph G(RA) where the set of nodes is 

multipliers for the n-m+1 dimensional recovery problem for 
n-m+l>2. This approach has also worked well to reduce the 

number of non-smooth iterations during recovery. 65 am j mcs 

Further, for track maintenance, initial feasible solutions 

are generated as follows. When a new scan of information (a a~{{z£&i 1 )\&a, j„*o t jj^o and there exists ziH2 . . . incz such 



N={zjta=l, ... ^1; in=i, , . . , Nm} 
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, that and j^}. 

Note that the nodes corresponding to zero index have not 
been Included in the above defined graph, since two vari- 
ables that have only the zero index in common can be 
assigned independently. Connected components of the graph 5 
are then easily found by constructing a spanning forest via 
a depth first search. (A detailed algorithm can be found in the 
book by Aho. Hopcroft and Ullman cited above). 
Furthermore, this procedure is used at each level in the 
relaxation, i.e.. is applied to each assignment problem [3.1] 10 
for n=3, . . . M. 

The original relaxation problem is decomposed first All 
relaxed assignment problems can be decomposed a priori 
and all recovery problems can be decomposed only after 
they are formulated. Hence, in the n-to-(n-l) case, we have is 
n-2 relaxed problems that can all be decomposed initially, 
and the recovery problems that are not decomposed (since 
they are all two dimensional). In the n-to-2 case, we have 
only one relaxed problem that can be decomposed initially. 
This case yields n-3 recovery problems, which can be 20 
decomposed only after they are formulated. 

The ever-increasing demand for sensor surveillance sys- 
tems is to accurately track and identify objects in real-time, 
even for dense target scenarios and in regions of high track 
contention. The use of multiple sensors, through more varied 25 
information, has the potential to greatly improve target state 
estimation and identification. However, to take full advan- 
tage of the available data, a multiple frame data association 
approach is needed. The most popular such approach is an 
enumerative technique called multiple hypothesis tracking 30 
(MHT). As an enumerative technique. MHT can be too 
computationally intensive for real-time needs because this 
problem is NP-hard. A promising approach is to utilize the 
recently developed Lagrangian relaxation algorithms for the 
multidimensional assignment problem; however, there are 35 
many other potentially good approaches to these assignment 
problems such as LP relaxation combined with an interior 
point method, GRASP, and parallelization. 

These data association problems are fundamentally 
important in tracking. Thus, our first objective is to present 40 
an overview of the tracking problem and the context within 
which these data association problems fit Following a brief 
review of the probabilistic framework for data association, 
the second objective is to formulate two models for track 
initiation and maintenance as multidimensional assignment 45 
problems. This formulation is based on a window moving 
over the frames of data. The first and simpler model uses the 
same window length for track initiation and maintenance, 
while the second model uses different lengths. As the 
window moves over the frames of data, one creates a 50 
sequence of these multidimensional assignment problems 
and there is an overlap region of frames common to adjacent 
windows. Thus, given the solution of one problem, one can 
develop a warm start for the solution to the next problem in 
the sequence, as shown hereinafter. Such information is 55 
critically important to the design of real-time algorithms. 

Overview of the Tracking Problem 

Tracking and data fusion are complex subjects of special- 
ization that can only be briefly summarized as they are 
related to the subject of this paper. The processing of track 60 
multiple targets using data from one or more sensors is 
typically partitioned into two stages or major functional 
blocks, namely, signal processing and data process. The first 
stage is the signal processing that takes raw sensor data and 
outputs reports. Reports are sometimes called observations. 65 
threshold exceedances. plots, hits, or returns, depending on 
the type of sensor. The true source of each report is usually 
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unknown and can be due to a target of interest, a false signal, 
or persistent background objects that can be moving in the 
field of view of the sensor. 

Two principal functions of data processing are tracking 
and discrimination or target identification. The cascrimina- 
tion or target identification function distinguishes between 
tracks that are for targets of interest and other tracks such as 
those due to background. Also, if there is enough 
information, each track is classified as to the type of target 
it appears to represent. Target identification and discrimina- 
tion will not be discussed further in this paper except to 
comment here on the use of attributes (including identifica- 
tion information) in the data association. 

There are' two types of attributes or features, namely, 
discrete valued variables and continuous valued variables. 
Available attributes of either or both types should be used in 
computing the log likelihoods or scores for data association. 
In discussing tracking in this paper, the attributes will not be 
mentioned explicitly but it should be understood that some 
reports and some tracks may include attributes information 
that should be and can be used in the track processing (both 
data association and filtering) if it is useful. 

Tracking Functions 

The tracldng function accepts reports for each frame of 
data and constructs or maintains a target state estimate, the 
variance-covariance matrix of the state estimation error, and 
refined estimates (or probabilities) of the target attributes. 
The state estimate typically includes estimates of target 
position and velocity in three dimensions and possibly also 
acceleration and other variables of interest. 

The tracldng function is typically viewed in two stages, 
namely, data association and filtering. Also, the process of 
constructing new tracks, called track initiation, is different 
from the process of updating existing tracks, called track 
maintenance. 

In track maintenance, the data association function 
decides how to relate the reports from the current frame of 
data to the previously computed tracks. In one approach, at 
most one report is assigned to each track, and in other 
approaches, weights are assigned to the pairings of reports 
to a track. After the data association, the filter function 
updates each target state estimate using the one or more 
(with weights) reports that were determined by the data 
association function. A filter commonly used for multiple 
target tracldng and data fusion is the well know Kalman 
filter (or the extended Kalman filter in the case of nonlinear 
problems) or a simplified version of it 

In track initiation, typically a sequence of reports is 
selected (one from each of a few frames of data) to construct 
a new track. In track initiation, the filtering function con- 
structs a target state estimate and related information based 
on the selected sequence of reports. The new track is later 
updated by the track maintenance processing of the subse- 
quent frames of reports. 

In some trackers, there is also a tentative tracking function 
for processing recently established tracks until there is 
enough confidence to include them in track maintenance. 
While for simplicity of presentation* tentative tracking will 
not be included in the processing discussed in this paper, the 
techniques discussed could readily include one or more 
tentative tracking functions. 

There have been numerous approaches developed to 
perform the data association function. Since optimal data 
association is far too complex to implement, good but 
practical sub-optimal approaches are pursued. Data associa- 
tion approaches can be classified in a number of ways. On 
way to classify data association approaches is based on the 
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number of data frames used in the association process. In reports is mot the same for these two functions. This is a 

single frame data association for track maintenance, "hard novel approach that is introduced in this paper, 

decisions" are made each frame as to how the reports are to Multiple Sensor Processing 

be related to the tracks. Some single frame approaches There are many advantages and many ways to combining 
include: individual nearest neighbor, PDA, JPDA. and gto- 5 data from multiple sensors. There are also many ways of 
bal nearest neighbor (sequential most probable hypothesis categorizing the different algorithmic architectures for pro- 
tracking), which uses a two dimensional assignment algo- cessing data from multiple sensors. One approach outlines 
rithm. In most single frame data association approaches, four generic algorithmic architectures. Two of these generic 
only one track per object is carried forward to be used in architectures are especially pertinent to this paper and are 
processing the next frame of reports. io summarized briefly. 

Multiple frame data association is more complex and In the Centralized Fusion algorithmic architecture, reports 

frequently involves purposely carry forward more than one are combined from the various sensors to form global tracks, 

track per target to be used in processing the next frame of This algorithmic architecture is also called Central Level 

reports. By retaining more than one track per target, the Tracking, Centralized Algorithmic Architecture, or simply 

tracking performance is improved at the cost of increased 15 the Type IV algorithmic architecture. In track maintenance 

processing load. The best known multiple frame data asso- for example, if a single frame data association is used, then 

ciation approach is an enumerative technique called multiple for data association a frame of reports from one sensor is 

hypothesis tracking (MHT) which enumerates all the global processed with the latest global tracks; then the global tracks 

hypothesis with various pruning rules. are updated by the filter function. After the processing of this 

As shown hereinafter, the log likelihood of the probability 20 frame of reports is completed, a frame of the reports from 

of a global hypothesis can be decomposed into the sum of another (or the same) sensor is processed with these updated 

scores of each track that is contained in the hypothesis. As global tracks. This process is continued as new frames of 

a consequence, the most probable hypothesis can be iden- data become available to the system as a whole, 

tified with the aid of a non-enumerative algorithm for the In Centralized Fusion, using the multi-dimensional 

assignment so formulated. 25 assignment algorithm for the data association with multiple 

Interface Between Track Initiation And Track Mainte- sensors is similar to the processing of data from a single 

nance sensor. Instead of using multiple frames of reports from a 

There are a number of ways that track initiation and track single sensor, the multiple frames of reports come from 

maintenance functions can interact. Two ways that are multiple sensors. The frames of reports from all the sensors 

especially pertinent to the methods of this paper will be 30 are ordered based on me nominal time each frame of reports 

discussed in this section. is obtained and independent of which sensor provided the 

One approach is to treat these two functions sequentially, reports. In this way the approaches discussed in this paper 

by first assigning reports to tracks and then using the can be extended to process reports from multiple sensors 

information that remains to form new tracks. A better using the Centralized Fusion algorithmic architecture, 
approach is to integrate both processes where initiating 35 The second pertinent algorithmic architecture is Track 

tracks compete for the new frame of reports on an equal Fusion. This approach is also called the Hierarchical or 

basis, i.e., at the same time. Federated algorithmic architecture, sensor level tracking or 

In the integrated approach, the data association for track simply the Type II algorithmic architecture. In track main- 
initiation and track maintenance processing are combined tenance for example, a processor for each sensor computes 
and conducted simultaneously. One assignment array is 40 single sensor tracks; these tracks are then forwarded to the 
created that includes the track scores for potential tracks for global tracker to compute global tracks based on data from 
both track initiation and track maintenance. The first dimen- all sensors. 

sion of the assignment problem includes only all the tracks After the first time a sensor processor forwards tracks to 
either created or updated in the processing of the prior the global tracker, then subsequent tracks for the same 
frames of reports. Each of the remaining dimensions accom- 45 targets are cross-correlated with the existing global tracks, 
modates one frame of reports. This track-to-track cross-correlation is due to the common 
After the assignment algorithm finds a solution, each of history of the current sensor tracks and the tracks from the 
the previously established tracks are updated with the report same sensor that were forwarded earlier to the global tracker, 
assigned to it from the second dimension of the assignment The processing must take this cross-correlation into account 
array. The remaining reports in the second dimension that 50 and there are a number of ways of compensating for this- 
were in the assignment solution are firmly established as cross-correlations. One method for dealing with this cross- 
track initiators. The unassigned reports in the second dimen- correlation is to decorrelate the sensor tracks that are sent to 
sion of the assignment array are discarded as false signals. the global tracker. There are a variety of ways to achieve this 
The processing of the next frame of reports repeats this decollation and some are summarized in Drummond, 
process using these updated tracks and the newly identified 55 1996, Signal Data Proc. of Small Targets 1995, 
initiators in the first dimension of the cost array for process- -2561:369-383. 

ing the new frame of reports. Details of this approach are Once me sensor tracks are decorrelated they can be 

discussed hereinafter. processed by the global tracker in almost the same way as 

The integrated approach just discussed, uses the same reports are processed. In the case of track fusion the asso- 
number of frames of reports for both track initiation and 60 ciation process is referred to as track (or track-to-track) 

track maintenance. The goal in using the multi-dimensional association rather than data (or report-to-track) association, 

assignment algorithm is to provide improved performance If the sensor tracks are decorrelated, the global tracker can 

while minimizing the amount of processing required. process the tracks from the various sensor processors in 

Typically, track initiation will benefit from more frames of much the same way that the global tracker of Centralized 
reports than will track maintenance. Thus a second approach 65 Fusion processes reports. Accordingly the methods 

that integrates track initiation and track maintenance is described in this paper can be readily extended to processing 

discussed hereinafter, wherein the number of frames of data from multiple sensors using either the Centralized 
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Fusion or Track Fusion or even a hybrid combination of both 



algorithmic architectures. 



Formulation of the Data Association Problem 
The goal of this section is to briefly outline the probabi- 
listic framework for the data association problems presented 
in this work. The technical details are presented elsewhere. 
Hie data association problems for multisensor and multitar- 
get tracking considered in this work are generally posed as 
that of maximizing the posterior probability of the surveil- 
lance region (given the data) according to 



Maximize {P(r=rfZ?+ l )^r*\ 



(3.1) 



where Z represents N+l data sets, y indices of the data 
(and thus induces a partition of the data). 

ns r 

TS yfi 



{z^ 1 z^** 1 } where each i k and z ik k can assume the 

values of 0 and Zq\ respectively. A partition of the data will 
refer to a collection of tracks of data wherein each report 
occurs exactly once in one of the tracks of data and such that 
all data is used up; the occurrence of dummy report is 
unrestricted The reference partition is that in which all 
reports are declared to be false. 

Next, under appropriate independence assumptions one 
can show 



/xr = r \z N+1 ) 



(33) 



L, . . . is a likelihood ratio containing probabilities for 
detection* maneuvers, and termination as well as probability 
density functions for report errors, track initiation and ter- 
mination. Define 



(3.4) 
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-continued 

1 if (z ii ..- (JV+J ) are assigned to as a track, 
0 otherwise. 



Then, recognizing that 



10 



fX)#\Z N +L) 



the problem (3.1) can be expressed as the N+l-dimensional 
assignment problem 



15 



20 



Minimize 



■I 



(3.5) 



Subject 



to y ... ^ 



a n = l. 



.Mi, 



partition, and P(TyiZ N+x ) is the posterior probability of a 
partition y N+l is defined below; however, this framework 
is currently sufficiently general to cover set packings and 25 
coverings. 

Consider N+l data sets 2(k)(k=l N+i). each 

consisting of M k actual reports and a dummy report Zq*. and 
let denote the cumulative data set defined by 

30 

2C*H^ ™ d ^=W\ - ■ ■ , ZCAW)}, (3.2) 

respectively. (The dummy report Zq* serves several purposes 
in the representation of missing data, false reports, initiating 
tracks, and terminating tracks.) In multisensor data fusion 
and multitarget tracking the data sets Z(k) may represent 
different classes of objects, and each data set can arise from 
different sensors. For track initiation the objects are reports 
that must be partitioned into tracks and false alarms. In our 
formulation of track maintenance, which uses a moving 
window, one data set will be tracks and remaining data sets 
will be reports which are assigned to existing tracks, as false 
reports, or to initiating tracks. 

We specialize the problem to the case of set partitioning 
defined in the following way. Define a *1rack of data" as 



for i k = 1 A/* and k * 2 N, 



where Cq 0 is arbitrarily defined to be zero. Here, each 
group of sums in the constraints represents the fact that each 
non-dummy report occurs exactly once in a 'track of data.** 
One can modify this formulation to include multiassign- 
ments of one. some, or all the actual reports. The assignment 
problem (3.5) is changed accordingly. For example, if z t * is 
to be assigned no more than, exactly, or no less than n t k 
times, then the "=1 M in the constraint (3.5) is changed to "s, 
=, = n ijt V* respectively. In making these changes, one must 
pay careful attention to the independence assumptions, 
which need not be valid in many applications. Expressions 
for the likelihood ratios can be found in the work of Poore, 
1994, Computational Optimization & AppL, 3:27-57, and 
the references therein. 
Track Initiation and Maintenance 
In this section we explain two multiframe assignment 
formulations to the track initiation and maintenance prob- 
50 lem. The continued use of all prior information is compu- 
tationally intensive for tracking, so that a window sliding 
over the frames of reports is used as the framework for track 
maintenance and track initiation within the window. The 
objectives are to describe an improved version of a simple 
55 method and then to put this into a more general framework 
in which track initiation and maintenance have different 
length moving windows. 
The Fhst Approach to Track Maintenance and Initiation 
The first method as explained in this section is an 
60 improved version of our first track maintenance scheme and 
uses the same window length for track initiation and main- 
tenance after the initialization step. The process is to start 
with a window of length N+l anchored at frame one. In the 
first step one there is only track initiation in that we assume 
65 no prior existing tracks. In the second and all subsequent 
frames, mere is a window of length N anchored at frame k 
plus a collection of tracks up to frame k. This window is 
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denoted by {k; k+1 k+N}. The following explanation 

of the steps is much like mathematical induction in that we 
explain the first step and then step k to step k+1. 
Track Maintenance and Initiation: Step 1. Let 



(4.1) 



be an enumeration of all those zero-one variables in the 
solution of the assignment problem (3.5) (i.e., z lj<2 t/m==1 ) 
excluding all the false reports in the solution (i.e., all those I0 
zero-one variables with exactly one nonzero index) and 
zero-one variables in the solution for which (i A ,i 2 >=(0.0). 
(The latter can correspond to tracks that initiate on frames 
three and higher.) These denote our initial tracks. 

Consider only the first two index sets in this enumeration 15 
(4.1) and add the zero index 1 2 =0 with the corresponding 
values of i, and i 2 being zero. Thus, the enumeration is now 
<i 1 (l 2 )J 2 (l 2 )}/ 2 ^0. The notation T^l^^^^) will be 
used for this pairing. Suppose now that the next data set, i.e., 
the (N+2)' A set, is added to the problem. 2 o 

To explain the costs for the new problem, one starts with 
the hypothesis that a partition yer* being true is now 
conditioned on the truth of the pairings on the first two 
frames being correct. Corresponding to the sequence {T 2 
(l 2 ),z l3 , . . . -Z tnf2 \* the likelihood function is then given by 2 $ 

L r 2 m = 



Next, define the cost and the corresponding zero-one vari- 
able 



2*2*3 +2 : 



1 if ^ + + 2 \ »s assigned to Track 7 2 (/ 2 ), 

0 otherwise 



respectively.. Then the track maintenance problem Maximize 
{Lyryer* assignment 

Minimize >>"' ^ <Vi»-'w+2*2fj~*42 



Subject to: 



> ... X ^'3-^2 = 1 * h = Uk Li, 

for i p = 1 M p said p = 4 tf + 2-1, 

**//+2-J 
'jV+2-I«° 



50 

-continued 

Track Maintence and Initiation: Step k. At the beginning 
of the k tb step, we solve the following (N+l)-dimensional 
assignment problem 



Minimize: ) c V*+i~ 'n 



(4.4) 



k= Q '1+1= 



Subject to: \ *** X 



= 1, 



for i p = l and /> = Jt + 2, ... , TV + k - 1, 

/*=<) i t+ i=0 



35 where for the first step 1, and Lj are replaced by \ i and 
respectively. The first index 1* in the subscripts corresponds 
to the sequence of tracks {T*(l fc )},/*=0 where T k (\ k )-{ z (/ l 

(l k ) z t A (ijt)} is a track of data from the solution of the 

problem prior to the formulation of (4.4) or if k=l then it is 

40 just the first data set in frame one. 

Suppose problem (4.4) has been solved and let the 
solution, i.e., those zero-one variables equal to one, be 
enumerated by 

{(WWiV**t(^i) U^.))}^- 1 < 4 * 5 > 

45 with the following exceptions. 

a. All zero one variables for which (l k d k+i )H0,0) are 
excluded. Thus, tracks that initiate on frames after the 
(k+1)"' are not included in the list. 

b. All zero-one variables whose subscripts have the form 
50 ijpO and exactiy one nonzero index in the remaining indices 

. . . 4*+*} ^ excluded. These correspond to false 
reports on frames p=k+l k+N. 

c. All variables (z v . . . ) for which (i* +J 

W)=(0,0 0) and Ul*M) in (4.5). In other words, the 

55 reports on the last N+l frames in the {T^IJ,^* 1 } are all 
dummy. A&y solution with this feature corresponds to a 
terminated track. 

Given the enumeration (4.5), one now fixes the assignments 
on the first two index sets in the list (4.5). The zero index 
60 1^=0 is added to the enumeration to specify (ijfc(0)a* +l (Q))= 
(0,0) that is used to represent false reports and tracks that 
initiate on frame k+2 or later, so that the enumeration (4.5) 
is now 
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(4.6) 



Then, for 1* +1 =1, . . . J^i* the 1, ** such track is denoted 
by {T w (WiWA(WWc+lJWi)}^ the (N+l)- 
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tuple {T^itW-C* WU1 den ° te ^f 

T k+I a^i) Plus a set of reports {z t J " ,z w h 

actual or dummy, that are feasible with the track T M (lfcnX 

The (N+l)-tuple {T M (0).z lk ^ 2 z w * ^> wiU 

denote a track that initiates in the sliding window. i.e.. on 
subsequent frames. A false report in the sliding window is 

one with all but one non-zero index ip for some p=k+2 

k+l+N in the (N+l)-tuple {T w (0)je^ 2 . - . . ^J +L+ % 
The corresponding hypothesis about a partition yeP 
being true is now conditioned on the truth of the L i+1 tracks 
existing at the beginning of the N-frame window. (Thus the 
assignments prior to this sliding window are fixed,) The 
likelihood function is given by 



10 



Given this notation for the tracks and partition of the data in 
the frames {k-L . . . . . . ±+J}. L r ^ will denote the 

accumulated likelihood ratio up to and including frame 
p(p^k) for a track that is declared as existing on frame k as 
a solution of the assignment problem. In this notation, the 
likelihood for T pk {l k ) and that of the association of 



} with track T^(l Jt ) is given by 



(4.7 a) 15 



respectively. 
20 The cost for the assignment of {z 



Next, define the cost and the zero-one variable by 

said at least one inert gas has an amount, on a molar basis, 
that is greater than the amount of said oxygen in said 25 
pressurized medium; 



... ■•■•0*> to 

track TjfcCljfe) and the corresponding zero one variable are 
given by 



= -ln[ L Tp j. , ijt , L t p+1 Uk) . , t (i t ), t+1 „. w ] 



(4.11) 



C 'i+l'i+2-'i+l+* 



(4.7b) 
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0 otherwise 



i if (C 



i 

0 otherwise 



respectively. Likewise, for costs associated with the false 
35 reports on the frames k-I to k and as associated with 



respectively, so that the track extension problem, which was 
originally formulated as Maximize {Lylyer*}, can be 
expressed as exactly the same multi-dimensional assignment 
in (3.4) but with k replaced by k+1. Thus, we do not repeat 
it here. 

The Second Approach to Track Maintenance and Initia- 
tion 

In the first approach the window lengths were the same f or 
both track maintenance and initiation. This can be inefficient 
in mat one can usually use a shorter window for track 
maintenance than for track initiation. This section addresses 
such a formulation. 

The General Step k. To formulate a problem for track 
initiation and maintenance we consider a moving window 

centered as frame k of length I+J+l denoted by [k-L 

k, . k+J]. In this representation, the window length for 

track maintenance is J and that for initiation. I+J+l. The 
objective will be to explain the situation at this center and 
then the move to the same length window at center k+1 

denoted by [k+l-L . . . Jt+1 Jc+l+J], i.e., by moving the 

frame one frame at time to the right. The explanation from 
the first step follows hereinafter. 

The notation for a track of data is 



able are given by: 



-z, M ] and me corresponding zero-one vari- 
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J' 

'i-r-'t'W'W 



1 * tV/ Vi+i Vj 1 

are assigned as a track, 
0 otherwise. 



(4.12) 
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For the window of frames in the range [k-L ... Jc], the 
easiest explanation of the partition of the reports is based on 
the definitions 

f i p | £ belongs to one of the tracks listed on frame A, 

te., to one of the T k {l t ) for some h - 1 U 

F^-tfl A/,]\?,jfc)U{0), 



(4.8) 



so that all of those reports not used in the tracks T k (lj) on 
6o frame p aire put into the set of available reports for the 
formation of tracks over the entire window. Thus, we 
formulate the assignment problem as 



where the index 1* is used for an enumeration of those 
reports paired together. We also use the notation T^l^) to 65 
denote the sequence of reports belonging to track T^y but 
restricted to frames prior to and including p. Thus, 



Minimize 



M M r 
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= l,/ 2 =l,/ 2 ^2. 



*2=0 i4=0 



1^ Mi 



= l,i 3 =1 W 3f 



zz-zz, 

/ 2 =0 J3=0 i p -t =(if p+l =i 



«/V+2 
W+2=° 



10 



for e jfc = I Mi and it = 2 N, 

z, r . Iaw e (0, 1) for all \ u ...Jn+i 



Suppose problem (4-14) has been solved and let the 
solution, i.e., those zero-one variables equal to one, be 



/ T enumerated by \ 

K/t + i> J *+« t, * +i> '"'* + ' < '* +i) '/ 4+l -i 



(4.15) 



with the following exceptions. 

a. All zero one variables in the second list for which 

(i k _,< . . . J A +1 M0 a0) ^ excluded - Tnus - tracks that 

initiate on frames after the (k+l)' A are not included in the list. 

b. All false reports axe excluded, i.e., all zero-one vari- 
ables in the second list whose subscripts have exactly one 
nonzero index. 

c. All variables z, 4tA+i ik J for which (i^, w)= 

(0,0 0) and z(p)nT Jt (i jt )=0 for p=k-K, . . . ,k where K^0 

is user specified. Thus the track T t (l Jfe ) is terminated if it is 
not observed over K+J+l frames. 

Given the enumeration (4.15), one now fixes the assign- 
ments on the all index sets up to and including the (k+1)"* 
index sets. 



decisions corresponding to the reports in our list of tracks 
prior to and including frame k+l=I+2. This defines the k for 
problem (4.4) and one can then continue the development by 
adding a frame to the window as in the general case. 

If I+J>N, then one possibility is to start the process with 
N+l frames, and assuming J^N, proceed as before replac- 
ing I by N-J for the moment, and continue to add frames 
without lopping off the first frame in the window until 
reaches a window of length I+J-f 1. Then we proceed as in the 
previous paragraph. 

If I+J<N, then one can solve the track initiation problem 
15 (3.5), formulate the problem with the center of the window 
at k+l=N+l-J, enumerate the solutions as above, and lop off 
the first N-J-I frames. Then, we proceed just as in the case 
I+J=N. 

20 A primary objective in this work has been to demonstrate 
how multidimensional assignment problems arise in the 
tracking environment. The problem of track initiation and 
maintenance has been formulated within the framework of a 
moving window over the frames of data. The solution of 
25 these NP-hard, noisy, large scale, and sparse problems to the 
noise level in the problem is fundamental to superior track 
estimation and identification. Thus, one must utilize the 
special structure in the problems as well take advantage of 
special information that is available. Since these moving 
30 windows are overlapping, there are some algorithm efficien- 
cies that can be identified and that take advantages of the 
overlap in the windows from one frame of reports to the 
next. Here is an example of the use of a primal solution of 
one problem to warm start the solution of the next problem 



35 



40 



in the sequence. 

Suppose we have solved problem (4.4) and have enumer- 
ated all those zero-one variables in the solution of (4.4) as 
in (4.5). Add the zero index 1^=0, so that the enumeration 
is 



for /jt + i = 1 I*, 

W^tfuiX ...,<<4+iX 

for /jt +1 = 1, Ljt+i + 1, .... 



(4,16) 45 



50 



With this enumeration one can define the cost by 
and the tw o dimensional assignment problem 



(5.1) 



(5.2) 



55 



Thus one can now formulate the assignment problem for 
the next problem exactly as in (4-14) but with k replaced by 
k+1. Thus, we do not repeat it here. 

The Initial Step. Here is one explanation for the initial 
step. First, assume that N=I+J. In this case, we start the track 
initiation with a solution of (3.5). Let 

\mMh) • • • ww^m**** 

be an enumeration of the solution set of (3.5), i.e., those 60 
zero-one variables z^m^ . . . w^ 1, including Zqo . _ 0 =1 
corresponding to 1 2 =0, but excluding all those zero-one 
variables that are assigned to one and correspond to false 
reports (Le., there is exactly one nonzero index in the 
subscript of z t%h all those zero-one variables that are 65 

assigned to one anefcorrespond to tracks that initiate on 
frames higher than 1+2. Then we fix the data association 



4>2 = Minimize £ £ A+0m+n ^'ui'w+tf s 

JWjt+l+// 

Subject to £ i.W =Ut+1 = 1 



Let w be an optimal or feasible solution to this two- 
dimensional assignment problem and define 
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Q 5 



55 



(5.4) 



1 if (i't+i *jt+Af) = <**+i (4+i )) 

"Wt+i+w = 1 for 4+i ~ l. — , i*+i 
or if (Wi +i+ *) = 0,0 
0 otherwise . 
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This need not satisfy the constraints in that there are usually 
many objects left unassigned. Thus, one can complete the 
assignment by using the zero-one variables in (4,4) with k 
replaced by k+1 with exactly one nonzero index correspond- 
ing to any unassigned object or data report 

For the dual solutions, the multipliers arising from the 
solution of the two dimensional assignment problem (5.3) 
corresponding to the second variable. Le„ {\^ k * L+ *}i k+H 
^k+i+K-Q These are good initial values for use in a relaxation 
scheme [11.12]. Finally, note that one can also develop a 
warm start for problem (4.14) in a similar fashion. 

Multidimensional assignment problems govern the cen- 
tral problem of data association in multisensor and multi- 
target tracking, i.e.. the problem of partitioning observations 
from multiple scans of the surveillance region and from 25 
either single or multiple sensors into tracks and false alarms. 
This fundamental problem can be stated as 



Minimize 



Subject to: > ... ^ ^...^ = 1, ,,-ijuj,, 
^ »a-=o 



for i p **\,...,M p and p = 2, ... , N - 1, 



zi v .., M e iO, 1} for all i u ... , i N , 

for = l M t aod £ = 2, ... , 

Zi y ..lfit+i e }0, 1) for all tj 



15 



20 



(1.1) 
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methods for solving this NP-hard problem optimally are 
enumerative in nature, with branch-and-bound being the 
most efficient; however, such methods are much too slow for 
real-time applications. Thus one must resort to suboptiraal 
approaches. Ideally, any such algorithm should solve the 
problem to within the noise level, assuming, of course, that 
one can measure this noise level in the physical problem and 
the mathematical method provides a way to decide if the 
criterion has been met 

There are many algorithms that can be used to construct 
suboptimal solutions to NP-hard combinatorial optimization 
problems. These include greedy (and its many variants), 
relaxation, simulated annealing, tabu search, genetic 
algorithms, and neural netork algorithms. For the three 
dimensional assignment problem, Pierskalla developed the 
tri-substitution method, which is a variant of the simplex 
method. Frieze and Yadegar introduced a method based on 
Lagrangian relaxation in which a feasible solution is recov- 
ered using information provided by the relaxed solution. Our 
choice of approaches is strongly influenced by the need to 
balance real-time performance and solution quality. 
Lagrangian relaxation based methods have been used suc- 
cessfully in prior tracking applications. An advantage of 
these methods is that they provide both an upper and lower 
bound on the optimal solution, which can then be used to 
measure the solution quality. These works extend the 
method of Frieze and Yadegar to the mulUdimensional case. 

Probabilistic Framework for Data Association. (ABP) 

The goaJ of this section is to explain the formulation of the 
data association problems that governs large classes of data 
association problems in centralized or hybrid centralized- 
sensor level multisensor/multitarget tracking. The presenta- 
tion is brief; technical details are presented for both track 
initiation and maintenance in the work of Poore. The for- 
mulation is of sufficient generality to cover the MHT work 
of Reid. Biackman and Stein and modifications by Kurien to 
include maneuvering targets. As suggested by Biackman. 
this formulation can also be modified to include target 
features (e.g.. size and type) into the scoring function. The 
recent work of Poore and Drummond significantly extends 
the work of this section to new approaches for multiple 
sensor centralized tracking. Future work will involve exten- 
sions to track-to-track correlation. 

The data association problems for multisensor and mul- 
titarget tracking considered in this work are generally posed 
as that of maximizing the posterior probability of the sur- 
veillance region (given the data) according to 



Maximize {P(r=^f)\^r*} 



(2-D 



where c 0 . . 0 is arbitrarily defined to be zero and is included 50 
for notational convenience. One can modify this formulation 
to include multiassignments of one. some, or all of the actual 
reports. The zero index is used in representing missing data, 
false alarms, initiating and terminating tracks. In these 
problems, we assume that all zero-one variables z t/ i N with 55 
precisely one nonzero index are free to be assigned and that 
the corresponding cost coefficients are well-defined. (This is 
a valid assumption in the tracking environment) Although 
not required, these cost coefficients with exactly one nonzero 
index can be translated to zero by cost shifting without 60 
changing the optimal assignment. Finally, this formulation is 
of sufficient generality to include the symmetric problem 
and the asymmetric inequality problem. 

The data association problems in tracking that are formu- 
lated as (1.1) have several characteristics. They are normally 65 
sparse, the cost coefficients are noisy and the problem is 
NP-hard, but it must be solved in real-time. The only known 



where Z N represents N data sets, y of the data (and thus 
induces a partition of the data), T5 



TS 



P(r y\Z N ) is the posterior probability of a partition y true 
given the data Z* The term partition is defined below; 
however, this framework is currently sufficiently general to 
cover set packings and coverings. 

Consider N data sets Z(k) (k=I, . . . JsT) each of M k reports 
{^*h/*=l. ^ d let z " denote the cumulative data set 
defined by 
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signments of one, some, or ail the actual reports. The 
assignment problem (2.5) is changed accordingly. For 
example, if z, * is to be assigned no more than, exactly, or no 
less than n,* times, then the *=1" in the constraint (2.5) is 
5 changed to' ^, =, =n, A V respectively. Modifications for 
group tracking and multiresolution features of the surveil- 
lance region will be addressed in future work. In making 
these changes, one must pay careful attention to the inde- 
pendence assumptions that need not be valid in many 
applications. 

Expressions for the likelihood ratios ,, N , can be 

found in the work of Poore. These expressions include the 
developments of Reid, Kurien, and Stein and Biackman. 
What's more, they are easily modified to include target 
features and to account for different sensor types. In track 
initiation, the N data sets all represent reports from N 
sensors, possibly all the same. For track maintenance, we 
use a sliding window of N data sets and one data set 



respectively. In multisensor data fusion and multitarget 
tracking the data sets Z(k) may represent different classes of 
objects, and each data set can arise from different sensors. 
For track initiation the objects are measurements that must 
be partitioned into tracks and false alarms. In our formula- 
tion of track maintenance, which uses a moving window 
over time, one data set will be tracks and remaining data sets 
will be measurements which are assigned to existing tracks, 
as false measurements, or to initiating tracks. In sensor level 
tracking, the objects to be fused are tracks. In centralized 
fusion, the objects may all be measurements that represent 
targets or false reports, and the problem is to determine 
which measurements emanate from a common source. 

We specialize the problem to the case of set partitioning 
defined in the following way. First, for notational conve- 
nience in representing tracks, we add a zero index to each of 
the index sets in (2.2). a dummy report Zq* to each of the data 
sets Z(k) in (2.2), and define a "track of data" as (z ( /*0 
where i* and z* can now assume the values of 0 and Zq 

respectively. A partition of the data will refer to a collection 20 conta iaing established tracks. The formulation is the same as 



10 



15 



of tracks of data wherein each report occurs exactly once in 
one of the tracks of data and such that all data is used up; the 
occurrence of dummy report is unrestricted. The dummy 
report z Q k serves several purposes in the representation of 
missing data, false reports, initiating tracks, and terminating 25 
tracks. The reference partition is that in which all reports are 
declared to be false. 

Next under appropriate independence assumptions one 
can show 
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P(V=<f>\Z N ) 



(2.3) 



above except that the dimension of the assignment problem 
is now N+l.@@@ 

To explain the costs for the new problem, one starts with 
the hypothesis that a partition y e T* conditioned on the truth 
of the pairings on the first two frames being correct. Cor- 
responding to the sequence {T 2 (l 2 ), z h z^}, the 

likelihood function is then given by 

r>= ]~| *w* +2 wheie 

L T , == 1. 



hi i N is a likelihood ratio containing probabilities for 35 
detection, maneuvers, and termination as well as probability 

density functions for measurement errors, track initiation Next, define the cost and the corresponding zero-one vari- 
and termination. Then if c it ^InL,, v able 



-4 



t'i,-4v)«r 



<2.4) 40 



Using (2.3) and the zero-one variable z i( . iy/ =l if Oi 

i N )e 7 and 0 otherwise, one can then write the problem (2. 1) 
as the following N-dimensional assignment problem: 



1 if <4 O* 

Zi 2 iy. -w+2 = assigned to Track T 2 (/ 2 ), 

0 otherwise 



(4.2b) 



Minimize ^jT Ci 1 „.i A ,Zi l -.i w 

Subject to £ Ziv-tN = 1 < f| = *• "* ' Ml) 

£ = I M ^ 

iih-iN 

Tj *i~to =1 
'i-Vi'p+r"'" 

(i p = l M p and p=2....N-l) 

2 «i — Jjy = 1 <to * x * ••*» 
'iV"to-l 

zi v .. lfl €{0, 1} for all ii 



respectively. Then the track maintenance problem Maximize 
{Z5) {LTiyer* assignment 
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where Co 0 is arbitrarily defined to be zero. Here, each 
group of sums in the constraints represents the fact that 65 
each-non-dummy report occurs exactly once in a **track of 
data." One can modify this formulation to include multias- 



Minimize } / " Zj C V* ~ •*ti+ iZt 2h'"to+2 
t^t^o to + 2=° 

Subject to: \ ... V zi 2 i y ~i N + 2 = U h = 1. '2, **• » 
£0 to + 2=o 

L ^™to*2 = 1 - 13 - 1 M ** 

^0,^0 to + 2=0 

for i k = 1 M k and £ = 2, ... A*, 
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-continued 

JW, u 

Zi^.M+i € iO, 1] for all 



Track Maintence and Initiation: Step k. At the beginning 
of the k** step, we solve the following (N+l)-dimensional 
assignment problem. 



60 



Then, for 1^=1 L w . the 1*^* such track is denoted 

by T^G^MV^^ (N+l)-tuple 

flW^iK^ • ■ • *w ' denote 

Tn-iflU) plus a set of reports {z^ ~ .z w 

actual or dummy, that are feasible with the track T A+1 (1, 



The (N+l)-tuple {T* +1 (0), z, 



(4.4) 



'}. will 

denote a track that initiates in the sliding window, i.e.. on 
subsequent frames. A false report in the sliding window is 

one with all but one non-zero index i^ for some p=k+2 

10 k+l+N in the (N+l)-tuple {T^(0). z,J* 2 z w ** A "*r- 

The corresponding hypothesis about a partition y e T* 
being true is now conditioned on the truth of the L fc+1 tracks 
existing at the beginning of the N-frame window. (Thus the 
assignments prior to this sliding window are fixed.) The 
15 likelihood function is given by 



Subject to: 



I- 



- 1. 4 = l, 



20 



= 1, 



for *^ = 1 Af ^ and p = A + 2, ...,*,+*+ 1, 

^W" 1 *-" € t°' 1} for ^ Ik * ik+i ii+N 



(4.7 a) 



25 



Next, define the cost and the zero-one variable by 



(47b) 



0 otherwise, 
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where for the first step l x and L : are replaced by i x and 
respectively. The first index 1* in the subscripts corresponds 
to the sequence of tracks {T^y where T^l^z^ 1 

(y % i *(y } is a track of data from the solution of the 

problem prior to the formulation of (4.4) or if k^l then it is 
just the first data set in frame one. 

Suppose problem (4.4) has been solved and let the 
solution, i.e.. those zero-one variables equal to one* be 
enumerated by 

with the following exceptions. 

a. All zero one variables for which (lA+iMO.O) are 
excluded. Thus, tracks that initiate on frames after the 
(k+l) rt are not included in the list. 

b. All zero-one variables whose subscripts have the form 
1^=0 and exactly one nonzero index in the re mainin g indices 
{W • * ■ *W} 316 excluded These correspond to false 
reports on frames p=k+l, . . - .k+N. 

c. All variables z^ . . . for which (ik+1 

ik+NM0,0 0) and i^l^H) in (4.5). In other words, the 

reports on the last N+l frames in the {T^ a J"" 1 } m 311 
dummy. Any solution with this feature corresponds to a 
terminated track. 

Given the enumeration (4.5), one now fixes the assignments 
on the first two index sets in the list (4.5). The zero index 
1^=0 is added to the enumeration to specify (\ k (0)i kJhl (0))= 
(0,0) that is used to represent false reports and tracks that 
initiate on frame k+2 or later, so that the enumeration (4.5) 
is now 



respectively, so that the track extension problem, which was 
originally f ormulated as Maximize {L^iyerS as exactly the 
same multi-dimensional assignment in (3.4) but with k 
replaced by k+1. Thus, we do not repeat it here. 
35 The Second Approach to Track Maintenance and Initia- 
tion 

In the first approach the window lengths were the same for 
both track maintenance and initiation. TCiis can be inefficient 
in that one can usually use a shorter window for track 
40 maintenance than for track initiation. This section addresses 
such a formulation. 

The General Step k. To formulate a problem for track 
initiation aud maintenance we consider a moving window 
centered as frame k of length I+J+l denoted by [k-L 
45 k, . . . Jt+J]. In this representation, the window length for 
track maintenance is J and that for initiation. I+J+l. The 
objective will be to explain the situation at this center and 
then the move to the same length window at center k+1 
denoted by [k+l-L . . . Jc+1. . . . Jc+l+J]. i.e.. by moving the 
50 frame one frame at time to the right. The explanation from 
the first step follows hereinafter. 
The notation for a track of data is 



(4-8) 



55 



where the index 1* is used for an enumeration of those 
reports paired together. We also use the notation T p J\f) to 
denote the sequence of reports belonging to track T^J but 
restricted to frames prior to and including p. Thus, 

WWWv^ • * ■ *w* w W> ior P*k-1.(4-S>) 

65 Given this notation for the tracks and partition of the data in 
the frames {k-L . . . X . . . Jc+J}. %.*(!*) will denote the 
accumulated likelihood ratio up to and including frame 
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p(p^k) for a track that is declared as existing on frame k as 
a solution of the assignment problem. In this notation, the 
likelihood for T_*(I A ) and that of the association of 
{zj~ l \,„ } with trackT t (J t ) is given by 
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-continued 

k k+J M r 

Subject to: £ £ ^<-/™Vw »W S 1 



respectively. 

The cc 
track Tjk( 
given by 



The cost for the assignment of {z- k+i 



10 



track Tjk(l A ) and the corresponding zero one variable are 



and 



* I*?. 



0 otherwise. 



(4.11) is 



20 



for (,ef 4 jt\|0) and q - A - /, ... Jfc, 

Z Z Z Z for ^ sBl m « 

and $ = Jt + 1 Jt + J, 

Vr-'tfc+i-W € l °' 11 for dl 



Z E^i+i-'i+i : 

p=t+l< r =0 



1, i q = 1, - 



for $ = A + 7, 



is assigned to track T A (l fc ), respectively. Likewise, for costs 
associated with the false reports on the frames k-I to k and 



as associated with \ z 



. ,z, M ) and the correspond- 2 5 



ing zero-one variable are given by: 

1 {if * t -j V/ zii+J \ 

^t-r"* k l M"'* M ~ 316 assig 1 *^ as a tract 
0 otherwise. 



Suppose problem (4.14) has been solved and let the 
solution, i.e., those zero-one variables equal to one, be 
enumerated by 



(4.12) 
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(4.15) 



with the following exceptions. 



For the window of frames in the range [k-L ... ,k], the 
easiest explanation of the partition of the reports is based on 
the definitions 



35 



a. All zero one variables in the second list for which 



(4.13) 



belongs to one of the tracks listed oa frame k. 



i.e., to one of the 7^/*) for some l k - 1, , L k 
^=({1 J#/>)\T,jfc)U|0h 



40 



45 



so that all of those reports not used in the tracks T^fl*) on 
frame p are put into the set of available reports F pjc for the 
formation of tracks over the entire window. Thus, we 
formulate the assignment problem as 



(i w ^ijt+iMO 0,0) are excluded. Thus, tracks that 

initiate on frames after the (k+l) th are not included in the list 

b. All false reports are excluded, i.e., all zero-one vari- 
ables in the second list whose subscripts have exactly one 
nonzero index. 

c. All variables z /|IJM . . . , for which (i t+1 , . . . ,W=(0, 

0 0) and Z(p)r>T A (l>{0} for p=k-K, . . . X where K^0 

is user specified. Thus the track T k (\ k ) is terminated if it is 
not observed over K+J+l frames. 

Given the enumeration (4.15), one now fixes the assign- 
ments on the all index sets up to and including the (k+l)' A 
index sets. 



= 1. •■• . U+l 



(4.16) 



Thus one can now formulate the assignment problem for 
the next problem exactly as in (4.14) but with k replaced by 
* t+j M r 60 k+1. Thus, we do not repeat it here. 

Minimize Z Z Z Z c *k-i ~ *kh+i - «*+/*-/•• Vjt+r-w + The Initial Step. Here is one explanation for the initial 

p=t ~ l 'p**?* ir=0 step. First, assume that N=I+J. In this case, we start the track 

initiation with a solution of (3.5). Let 

^.imt+iw ^ aQ enymera tion of the solution set of (3.5), i.e M those 

zero-one variables z,^^ . . . ^(^=1, including 200 , 0=! 
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corresponding to i 2 =0. but excluding all those zero-one 
variables that are assigned to one and correspond to false 
reports (i.e., there is exactly one nonzero index in the 
subscript of z 



■W-t-I) 



all those zero-one variables that are 
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assigned to one and correspond to tracks that initiate on 
frames higher than 1+2. Then we fix the data association 
decisions corresponding to the reports in our list of tracks 
prior to and including frame k+l=I+2. This defines the k for 
problem (4.4) and one can then continue the development by 
adding a frame to the window as in the general case. 

If I+J>N. then one possibility is to start the process with 
N+l frames, and assuming J=!N. proceed as before replac- 
ing I by N-J for the moment, and continue to add frames 
without lopping off the first frame in the window until 15 
reaches a window of length I+J+L Then we proceed as in the 
previous paragraph. 

If I+J<N. then one can solve the track initiation problem 
(3.5), formulate the problem with the center of the window 
at k+l=N+l-J, enumerate the solutions as above, and lop off 
the first N-J-I frames. Then, we proceed just as in the case 
I+J=N. 

A primary objective in this work has been to demonstrate 
how multidimensional assignment problems arise in the 2 $ 
tracking environment. The problem of track initiation and 
maintenance has been formulated within the framework of a 
moving window over the frames of data. The solution of 
these NP-hard, noisy* large scale, and sparse problems to fee 
noise level in the problem is fundamental to superior track 30 
estimation and identification. Thus, one must utilize the 
special structure in the problems as well take advantage of 
special information that is available. Since these moving 
windows are overlapping, there are some algorithm efficien- 
cies that can been identified and that take advantages of the 35 
overlap in the windows from one frame of reports to the 
next Here is an example of the use of a primal solution of 
one problem to warm start the solution of the next problem 
in the sequence. 

Suppose we have solved problem (4.4) and have enumer- 
ated all those zero-one variables in the solution of (4.4) as 
in (4.5). Add the zero index 1 A+1 =0. so that the enumeration 
is 



With this enumeration one can define the cost by 
and the two dimensional assignment problem 

Subject to: £ 



Let w be an optimal or feasible solution to this 
dimensional assignment problem and define 
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(5.4) 



**** *Wi+J+JV ~ f ° r SOSOe = l * 1 

or if (^ + i,ii+i+/v) = (0, OX 
0 otherwise. 



This need not satisfy the constraints in that there are 
usually many objects left unassigned. Thus, one can com- 
plete the assignment by using the zero-one variables in (4.4) 
with k replaced by k+1 with exactly one nonzero index 
corresponding to any unassigned object or data report 

For the dual solutions, the multipliers arising from the 
solution of the two dimensional assignment problem (5.3) 
corresponding to the second variable, i.e.. \u k * i * N } 
/Wn-jQ. These are good initial values for use in a relaxation 
scheme [11.12]. Finally, note that one can also develop a 
warm start for problem (4.14) in a similar fashion. 

Multidimensional assignment problems govern the cen- 
tral problem of data association in multisensor and multi- 
target tracking, i.e.. the problem of partitioning observations 
from multiple scans of the surveillance region and from 
either single or multiple sensors into tracks and false alarms. 
This fundamental problem can be stated as 



Minimize 



2] -2 



(1.1) 



'1= 
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p = 1 M p and p = % .... AT-1, 



(5.1) 



(5.2) 50 



Zt y ..t S € {0, 1} for all fi i N , 

for i k = 1 Aft and U2,...,N, 

V-Vn €{0 ^ 11 foraU 'l 



where c 0 0 is arbitrarily defined to be zero and is included 
tS.3) f or notational convenience. One can modify this formulation 
55 to include multiassignments of one, some, or all of the actual 
reports. The zero index is used in representing missing data, 
false alarms, initiating arid terminating tracks. In these 
problems, we assume that all zero-one variables z tj i^ with 
precisely one nonzero index are free to be assigned and that 
the corresponding cost coefficients are well-defined, CThis is 
a valid assumption in the tracking environment.) Although 
not required, these cost coefficients with exactly one nonzero 
index can be translated to zero by cost sriifting without 
65 changing the optimal assignment Finally, this formulation is 
two- of sufficient generality to include the symmetric problem 
and the asymmetric inequality problem. 
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The data association problems in tracking that are formu- Consider N data sets Z(k) (k=I, . . . JST) each of reports 

lated as (1.1) have several characteristics. They are normally {z lk k \ tk Mk =[> and let Z N denote the cumulative data set 

sparse, the cost coefficients are noisy and the problem is defined by 
NP-hard, but it must be solved in real-time. The only known 

methods for solving this NP-hard problem optimally are 5 Z^H^*},"*^ andZ*={Z(i), . ■ - *Z(A0} (2.2) 
enumerative in nature, with branch-and-bound being the 

most efficient; however, such methods are much too slow for respectively. In multisensor data fusion and multitarget 
real-time applications. Thus one must resort to suboptimal tracking the data sets Z(k) may represent different classes of 
approaches. Ideally, any such algorithm should solve the objects, and each data set can arise from different sensors, 
problem to within the noise level, assuming, of course, that 10 p or initiation the objects are measurements that must 
one can measure this noise level in the physical problem and be partitioned into tracks and false alarms. In our formula- 
the mathematical method provides a way to decide if the tion 0 f maintenance, which uses a moving window 
criterion has been met. over time, one data set will be tracks and remaining data sets 
There are many algorithms that can be used to construct w ju be measurements which are assigned to existing tracks, 
suboptimal solutions to NP-hard combinatorial optimization 15 as f a i se measurements, or to initiating tracks. In sensor level 
problems. These include greedy (and its many variants), tracking, the objects to be fused are tracks. In centralized 
relaxation, simulated annealing, tabu search, genetic fusion, the objects may all be measurements that represent 
algorithms, and neural netork algorithms. For the three targets or false reports, and the problem is to determine 
dimensional assignment problem, Pierskalla developed the which measurements emanate from a common source, 
tri-substitution method, which is a variant of the simplex 20 We specialize the prob i em t0 tne case of set partitioning 
method. Frieze and Yadegar introduced a method based on defined in me following way. First, for notational conve- 
Lagrangian relaxation in which a feasible solution is recov- nience in re p res enting tracks, we add a zero index to each of 
ered using information provided by the relaxed solution. Our the index sets in (2 .2), a dummy report to each of the data 
choice of approaches is strongly influenced by the need to gets Z(k) in (2 2) and define a *^ ack of data » as ( z ") 
balance real-time performance and solution quality. 25 where @^ and z * can now assume the values of 0 and 
Lagrangian relaxation based methods have been used sue- ^ respective i y . A * partition of the data will refer to a 
cessfuUy in prior tracking applications. An advantage of collection of tracks of data wherein each report occurs 
these methods is that they provide both an upper and lower exact i y once in one of the tracks of data and such that all data 
bound on the optimal solution, which can then be used to i$ uged up; me occurrence G f dumm y report is unrestricted, 
measure the solution quaUty. These works extend the 30 The dummy report ^ serves severa i purposes in the rep- 
method of Frieze and Yadegar to the multidimensional case. reS entation of missing data, false reports, initiating tracks. 
Probabilistic framework for Data Association. (ABP) and terminating tracks. The reference partition is that in 
The goal of this section is to explain the formulation of the which all reports are declared to be false, 
data association problems that governs large classes of data Next under appropria te independence assumptions one 
association problems in centralized or hybrid centralized- " can show 
sensor level multisensor/multitarget tracking. The presenta- 
tion is brief; technical details are presented for both track ^ (2 3) 
initiation and maintenance in the work of Poore. The for- ~^ N =l y = f] V 
mulation is of sufficient generality to cover the MHT work Ar = y°(Z ) ^...i^er 
of Reid, Blackman and Stein and modifications by Kurien to 40 
include maneuvering targets. As suggested by Blackman, 

this formulation can also be modified to include target L if _ i N is a likelihood ratio containing probabilities for 
features (e.g., size and type) into the scoring function. The detection, maneuvers, and termination as well as probability 
recent work of Poore and Drummond significantly extends density functions for measurement errors, track initiation 
the work of this section to new approaches for multiple 45 and termination. Then if c tj i^^JnL^ . . . i*, 
sensor centralized tracking. Future work will involve exten- 
sions to track-to-track correlation. I P{y\Z N )\ y, < 2 - 4 > 

The data association problems for multisensor and mul- l Ay 0 \z N )\ ~ (Ji y N)f£y Cii '" iN ' 
titarget tracking considered in this work are generally posed 
as that of maximizing the posterior probability of the sur- 
veillance region (given the data) according to Using (2 j) and the zero-one variable =11^ N ) 

Maximize {/ , (T^y*Z ft ')^r*} (2.1) [2 5, 

€ y * l 

where Z N represents N data sets, y of tne ( and tnus 55 ^ 
induces a partition of the data), T5 



50 



Minimize £ ^i^n^i^n 



Y$ y ° 60 Subject to £ „j N ~ 1 (i t = 1, M t ) 



P(r tIZ^) is the posterior probability of a partition 7 true 

given the data The term partition is defined below; 65 2j = 1 

however, this framework is currently sufficiently general to 'i-ViVu 
cover set packings and coverings. 
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-continued 

(«, = !,...,*/, ana p = 2,.../V-l) 



z 

'1*2 -*N-l 



10 



15 



Zi x ... w € {0, 11 for all *! i*- 



where Cq 0 is arbitrarily defined to be zero. Here, each 
group of sums in the constraints represents the fact that 
each-non-dummy report occurs exactly once in a "track of 
data." One can modify this formulation to include multias- 
signments of one, some, or all the actual reports. The 
assignment problem (2.5) is changed accordingly. For 
example, if z, 4 * is to be assigned no more than, exactly, or no 
less than n, * times, then the "=1" in the constraint (2.5) is 
changed to =. =, respectively. Modifications for 

group tracking and multiresolution features of the surveil- 
lance region will be addressed in future work. In making 20 
these changes, one must pay careful attention to the inde- 
pendence assumptions that need not be valid in many 
applications. 

Expressions for the likelihood ratios L h i^, can be 
found in the work of Poore. These expressions include the 
developments of Reid. Kurien. and Stein and Blackman. 
What's more, they are easily modified to include target 
features and to account for different sensor types. In track 
initiation, the N data sets all represent reports from N 
sensors, possibly all the same. For track maintenance, we 
use a sliding window of N data sets and one data set 
containing established tracks. The formulation is the same as 
above except that the dimension of the assignment problem 
is now N+l. 

Overview of the Lagrangian Relaxation Algorithm. 
(ABP) 

Having discussed the N-dimensional assignment problem 
(1.1), we now turn to a description of the Lagrangian 
relaxation algorithm. The algorithm will proceed iteratively 40 

for a loop k=l N-2. At the completion, there remains 

one two-dimensional assignment problem that provides the 
last step which yields an optimal (sometimes) or near- 
optimal solution to the original N-dimensional assignment 
problem. In step k of this loop (summarized in Section 4) 4 ? 
one starts with the following (N-k+l)-dimensional assign- 
ment problem with one change in notation. If k=l. then the 
index notation l x and L : are to be replaced by i L and M r 
respectively. 
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Subject to 

for 1, = 1 M p and P = k+Z .... N-l, 



(3.1) 



To ensure that a feasible solution of (3.1) always exists for 
is a sparse problem, all variables with exactly one nonzero 
index (i.e., variables of the form z 10 ^ 



L for L. 



andz^ 0 . 0 "~*+ l for y=l, . . . . and p=k+12. . . . 
,N) assumed free to be assigned with the corresponding cost 
coefficients being well-defined. This assumption is valid in 
the tracking environment. 

Subsection 3.1 presents some of the properties associated 
with the Lagrangian relaxation of (3.1) based on relaxing the 
last (N-k)-s;ets of constraints to a two-dimensional one. A 
new approach to the problem of recovering a high quality 
feasible solution of the original (N-k+1) -dimensional prob- 
lem given a feasible solution (optimal or suboptimal) of the 
relaxed two-dimensional problem is described hereinafter. A 
summary of the relaxation algorithm is given hereinafter, 
following which we establish the maximization of the 
Lagrangian dual (an important aspect of the relaxation 
procedure) to be an unconstrained nonsmooth optimization 
problem and then present a method for computing the 
subgradients. 

The Lagrangian Relaxed Assignment Problem. The 
(N-k+l)-dimensional problem (3.1) has (N-k+1) sets of 
constraints. A (M p +l)-dimensional multiplier vector (Le.. xP 
e R Mp 1) associated with the p-th constraint set will be 

denoted by vP=(\xgnf u^ p ) r with u/sO being fixed 

for each p=k+2 N and included for notational 

convenience only. Now, the (N-k+l)dimensional assignment 
problem (3.1) is relaxed to a two-dimensional assignment 
problem by incorporating (N-k-1) sets of constraints into 
the objective function via the Lagrangian. Although any 
(N-k-1) constraint sets can be relaxed, we choose the last 
(N-k-1) sets of constraints for convenience. The relaxed 
problem is 
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-.«")= Minimise (fa-fid"" 4 *'; «** 2 - •«") = 

M ~ Z <£ , -*<& , -«.* 



(3.2) 



ZZ< 

/>=*+2 < p =0 



Minimize 



L 
z z< 

Subject ,o X ^U-^U^l..-.^ 



One of the major steps in the algorithm is the maximization Then. 

25 4y,_ A+I «/** 2 , 



Of <I>* 



,u") with respect to the multipliers 



(u* +2 , . . . ,u*). It turns out that <!V* +1 is a concave, 
continuous, and piecewise affine function of the multipliers 
(u* +2 , . . . ,u"), so that the maximization of is a 

problem of nonsmooth optimization. Since many of these 
algorithms [???] require a function value and a subgradient 

of <J>yv-ik+i at anv required multiplier value (u A+2 u*'), we 

address this problem in the next subsection. We note, 
however, that there are other ways to maximize and 
the next subsection just addresses one such method. 
Properties Lagrangian Relaxed Assignment Problem 
For a function evaluation of ^^+1, we show that an 
optimal (or suboptimal) solution of this relaxed problem 
(3.2) can be constructed from that of a two-dimensional 
assignment problem. Then, the nonsmooth characteristics of 
^Ar-it+i are addressed, followed by a method for computing 
the function value and a subgradient. 

Evaluation of <J>awh-i* Define for each (l^i^i) an index 
0** - • •. Ja/KJ^OA-mX • . . jj/Wm)) and a new cost 



30 



35 



Minimize ^.^(z 2 ; u t+2 , ... , u N ) = 
Subj. To £ 4^ = 1,/* = !,...^ 



(3.4) 



.=0 



i t =a 

Ait+i e ^ l] for a11 



function c. 



'Vah 



by 



40 As an aside, two observations are in order. The first is that 
the search procedure needed for the computation of the 
relaxed cost coefficients in (3.3) is the most computationally 
intensive part of the entire relaxation algorithm. The second 
is that a feasible solution z N ~ k+i of a sparse problem (3.1) 

45 yields a feasible solution z 2 of (3.4) via the construction 



arg mLn< 



(3.3) 



I ,=i+2 



f 1, if # k J™i M ...i N =1 
( 0, otlicrwisc . 



for some {ik+2....j N % 



1 M ? and p = k + 2,..., N 



for (/i,ii + i)*(0,0) 



1 



50 



thus, there are generally solutions other than the one nonzero 
index solution. 

The following Theorem 3. 1 gives a method for evaluating 
55 **V-a+i a*"* states that an optimal solution of (3.2) can be 
computed from that of (3.4). Furthermore, if the solution of 
either of ithese two problems is e-optimal, then so is the 
other. The converse is contained in Theorem 3.2. 
Theorem 3.1. Let (o 2 dimensional assignment problem 
60 (3.4) and define co^ +l 



Given an index pair (l k d k+ i),Q k +^ . . . , j*) need not be 65 
unique, resulting in the potential generation of several 
subgradients (3.11). 



and </ t ,t*+i)*<0,( 
h*ti M ..Jt = 0 if f'jt+2, ... *n) * (it+2 Jn) 



(3.5) 
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-continued 

and (4>4+i)*(0,( 



= J 



Jt+1 

*k+2- J N ' 



if A' 



'jfc+2~W 



►I" 

/w*+2 



p=ktl 



>0. 



Then or d problem and <b N - k +xyd 
<*>amc-iY^ V k+2 x>"A addition, co, ar=* +1 O^jfyV^ 
O^iYD* 2 v N A 

Pro"ol Let to"-** 1 co 2 O^yco™ a/ +2 a>*A <D 
Y/^+i^*" 1 " 2 ^ objective function values of (3.2) and (3.4), 
respectively. Direct verification shows that co* * l in (3.2) 
and ^A^t+iTCO^ 1 \)* +2 v" 4>^ + i7C0 2 D*" 2 . . . .^remain- 
der of the proof, assume that CO" y A 
Letx*"*" 1 yA 

^o-l if (4,4+i) = (0,0) and 



for (4.4+1)* (0,0); 



(3.6) 



WV-A+1 



^ «^ <0 for some (4 +2 , ijv); 



^oo=0 if (4,4+i) = (0.0) and 



f + H "*p > 0 for ^ (i * +2 

;wfc+2 



Note mat x 2 Y A 



(or 



y A (0 (0 



optimal solution of (3.2), 



D A <JW^(co 2 -o^ 2 d A co Y A 

With the exception of one equality being converted to an 
inequality, the following theorem is a converse of Theorem 
3.1. 

Theorem 3.2. Let co*-** 1 y A and define co 2 
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The Nonsmooth Optimization Problem. Next, we address 
the nonsmooth properties of the function ^V^-nYU** 2 v N A 
explained in the following theorem 

Theorem 3.4 Let <LV^ y A ^WxY*"^ 1 ) be the object 
function value of the (N-k+l>dimensional assignment 
problem in equation (3.1) corresponding to a feasible solu- 
tion z""^ 1 of the constraints in (3.1). and let z^ 1 be an 
optimal solution of (3.1). Then. ^V^iY^ 2 v*A is piece- 



wise affine. concave and continuous in {v kJr * v 



. -t+1 

V*+i4+a-w 



for (4, 4+i)* (0,0); 



(3.7) 



ojgo = 1 if (4, 4+i) = (0, 0) and 

c «?!2™^ + Z ^ £ ° for some (i *+ a iAf); 

^Jfc+2 

4o=0 if (4*4+i) = (0,0) and 



45 



and 



10 



Most of the algorithms for nonsmooth optimization are 
based on generalized gradients, called subgradients. given 
by the following definition for a concave function. 

Definition 3.5. At u=(u* +2 u") the set 5<*V_* +I (u) is 

called a subdifferential of $ 



/•vV-i+l 



^ vF tp >0 for all (4+2, 



p*k+2 



AWfc+1 



and is defined for the 

concave function <5>yv_^i(u) by 

i(«)^ r (o>-«) for all coeR*" +J+1 * . . . jcR"*^ 

where are all permanently fixed. (Recall that 

these were used for notational conveience only.) A vector 
ge34>A/_jt +1 is called a subgradient. 

There is a large literature on such problems, and we have 
tried a variety of methods including subgradient methods 
and bundle methods. Of these, we have determined that for 
a fixed number of nonsmooth iterations (e.g.. twenty), the 
bundle trust method Schramm and Zowe provides excellent 
quality solutions with the fewest number of function and 
subgradient evaluations, and is therefore our currently rec- 
ommended method. 

An Algorithm for Computing and a Subgradient 

Most current software for maximizing the concave function 
requires the value of the function and a subgradient 

at a point (u* +2 u*). Based on the previous two sections, 

this can be summarized as follows. 

1. Starting with problem (3.1). form the relaxed problem 
(3.2); 

2. To solve (3.2) optimally, defined the two dimensional 
assignment problem (3.4) via the transformation (3.3); 

3. Solve the two-dimensional problem (3.4) optimally; 

4. Reconstruct the optimal solution, say of (3.2) 
via equation (3.5) as in Theorem 3.1; 

5. Compute the function 



<Wn<^ +2 «*) = 



2 

4 **+!♦— 1 n 



(3.10) 



- N-k+i 

■'AF^Vt+l- 
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N Up 
p=i+2 lp*Q 



Vfc+l-V-l'/H-i*"** 



- N-k+l 



"Then (0 2 Y A <D^ +1 y<o"-* +1 ; i>* +2 ^^jf^lf jf* 6. A subgradient is given by 
v A co"-** 1 y A co (3.4), *^,Yto" ^ v N A= & & J 

<*W<0 2 ^ »*A and (D^, (tf~ 2 v»A 4>^ +1 yi) A+2 _ 



Z. N-k+1 . 
^/4+r"^" 1 

44+i "^lV+i*" 1 * 



Proof: Let co v -* +1 co 2 * J ^ 1 (fli w ^ M ^ 2 v N A and t2 N 

Qft-k+ii® 2 vk+2 objective function values of (3.2) and 60 b^-^ii^ 

(3.4), respectively. Direct verification shows that CO 2 (3.4) bu f p 

and (SWiW^ 1 ^ 2 ^^-W® 2 ^ ^ A remain - 
derof me prcof. assume mat ©^^Y A and construct© 2 co for i P = i, m p and p = Jt+2,..., N. 

by replacing co^~* +1 y A co AWm " 1 co 

&N^ifX-*+W* 2 ^A=d^(co 2 o>* +2 x> N A=z 65 

^AMfc+iY 05 ^^ 1 * v ** 2 ^A CO^ inequality is in fact an Several remarks are in order. First, the optimal solution of 
equality. This proves the theorem, the nunimiation problem (3.2) is required before one can 
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remove the minimum sign, replace z*"** 1 by the minimizer 
^N-k+i an( j diff eren ti a te with respect to u,^ to obtain a 
subgradient as in (3.11). If fy A, ~ k+l is an approximate solu- 
tion of (3.2), then the subgradient and function values are 
only approximate with accruacy depending on that of 5 
^jv-A+i Although one can evaluate the sums in (3.10) and 
(3.11) in a straightforward manner, another method is based 
on the following observation. Given the multiplier vector 
«*+2 ^t {(Wl^OJ^itWi)}^ 1 ^. be an enu- 
meration of indices {1*4*+!} of w 2 (or the first 2 indices of 1Q 
w*~* +1 constructed in equation 3.5)) such that 
W «WWW^1 for (WW4 + A + iM0,0) and (lA+i)- 
ijN-i(Wi)MO'O) for regardless of whether Wq^I or 

not. (The latter can only improve the recovered feasible 
solution.) The evaluation of the bracketed quantity in (3.11) i5 

for a specific i p >=l and p=k+2 N is one minus the 

number of times the index value i p appears the (p-k+1)*' 
position of the (N-k-t-l)tuple in the list ({ljfc(ljfc +1 )4jt+i(U+i)' 
1** ■ * - o^W/^o. 

Finally, we have presented a method for computing one 20 
subgradient. If ^ Af ~^ 1 is the unique optimal solution of 
(3.2), then CV^jfu) is differentiable at u, g is just the 
gradient at u, and the subdifferential oO^+^u^jg} is a 
singleton. If, corresponding to u,vV /v ~* +i is not unique, then 
there are finitely many such solutions, say ^ 

{vy^V) ,\v"~ k+1 (m)} with respective subgradients 

{g(l) g(m)}, the subdifferential50(u)is the convex hull 

of {g(l), g(m)}. These nonunique solutions arise in two 

ways. First, if the optimal solution of the two dimensional 
assignment problem is not unique, then one can general all 30 
optimal solutions of the two dimensional assignment prob- 
lem (3.4). Corresponding to each of these solutions and to 
each index pair CW^m each solution, compute the indices 
* * * Av) i n (3*3). If these j*s are not unique, then we can 
enumerate all the possible optimal solutions iv* - *** 1 of (3.2). 35 
Given these solutions, one can then compute the correspond- 
ing subgradients from (3.11). 

In most nonsmooth optimization algorithms, one uses an 
epsilon-subdifferential 

Definition 3.x. At u=(u* +2 u") the set S e <I> // _ Jfc+1 (u) is ^ 

called an epsilon subdifferential of and is defined for 

the concave funtion ^> Ar _* +1 (u) by 

i(u)£g T (v>~u)te for all wR Mt3+, x . . . xR*"* 1 }. (3.9) 

45 

where gtf=tof*to(f=Q are all permanently fixed. (Recall that 
these were used for notational convenicence only.) A vector 
ged^^jt^u) is called a subgradient. 

HOW DO YOU COMPUTE gcBO^^u) 

The Recovery Procedure. The next objective is to explain 50 
a recovery procedure, i.e., given a feasible (optimal or 
suboptimal) solution to 2 y A y (a* - ** 1 y A via Theorem 3.1), 
generate a feasible solution z""* +1 of equation (3. 1) which 
is close to co 2 specified. We first assume that no variables in 
(3.1) are presassigned to zero; this assumption will removed 55 
shortly. The difficulty with the solution & N ~ k+l to it need not 
satisfy the last (N-k-1) sets of constraints in (3.1). (Note, 
however, that if co 2 (3.4) and a)*"*"*" 1 co relaxed constraints, 
then co*"* 4 " 1 y A A co recovery proceudre described here is 
designed to preserve the 0-1 character of the solution to 2 y 60 
A If fi> /Aiui 2 -l and 1^5*0 or i^^O, the corresponding feasible 
solution V*~* +l of (3.1) is construed so that ^ 0 
o N ~** l =l for some (i A+2 , . . . i^). By this reasoning, variables 
of the form can be assigned a value of one 

in the recovery problem only if Crtoo 2== *- However, variables 65 
Zoo^ . x jy A '~* +1 wm ^ be treated differently in the recovery 
procedure in that they can be assigned either zero or one 



independent of the value of a^ 2 . This increases the feasible 
set of the recovery problem, leading to a potentially better 
solution. 

Let UlitUjk+iW^iO^)}^!*** 1 ^ oe an enumeration of 
indices L +l \ of or y co*"** 1 equation (5)) such that 

for OA+iXWWMO.O) and 
i*+i(Wi3W0,0) for 1^=0 regardless of whether 00ocf = l or 
not. To define the N-k dimensional assignment problem that 
resores feasibility, let 



for k = 1 



(3.12) 



for k > 2 and i k = 0, L k 



where 

L . . w=*>Ui* • • * • • • 0W4*i)) - - - )X3-13) 

for m=2 k+1 and where o denotes function composition. 

For example l 2 ^2(h) and I234J2 0 WUM^WUX 

Then the (N-k)-dimensional assignment problem that 
restores feasibility is 



(3.14) 



Minimize £ 

Subject to £ <;^ +2 ... w = 1. l M - 1 

Z Ci'a* -in = = 1 A/t+2 - 

Z Oi+J-'Af " *• 

for /> = 1, ...,M p and p = JC + 3, N- I, 

Ctw"'* E {0 ' I( for aU **** ^ '** 



Let Y be an optimal or feasible solution to this N-k 
dimensional assignment problem. The recovered feasible 
solution z N is defined by 



N 



(3.14) 



[ 1 if /1 = k{l 2 ... mi. h = i 2 (h». t+i ), '3 * hih -k+i)* — 
't ~k{tkk+i\ k+i = U+iiU+i) and Xt i+l i k+2 - i N = t, 
\0 otherwise. 



Said in a different way, the recovered reasible solution 7? 

corresponding to the multiplier set \u k+2 u*} is then 

defined by 



•n ('2 - t+l y2<*2". A+l ^3%-.. t+i Vjt tfu+i V M U M U l+2 -iff : 



0 otherwise 



where l m k+l is defined in (3.13) and o denotes function 
composition. 



15 
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The next objective is to remove the asumption that all cost 
coefficients are defined and ail zero-one variables are free to ^(*2...jv-i)k<fc...ff-iH3!^« * t3J9) 

be assigned. We first note that the above construction of a f i if y / — 

reduced dimension assignment problem (3.14) is valid as 1 *-i ' 

long as all cost coefficients c"~ k are defined and all zero-one 5 1° otherwise, 

variables in are free to be assigned. Modifications are 

necessary for sparse problems. If the cost coefficient c ww whefe ^ ^ ^ defined in (3 D) and o denotes function 
'mCM*. • • - 1* +1 is_well-defined and the zero-one variable composition. To complete the description, let {ft-iOWa* 
^(W^KVx^- ^ x is free to be assigned to zero or one. 10 (i^)}, **=0be an enumeration of indices {W^} of Ysuch 
then define c Ww . ^^ +lWW , i+2 . ^ as in ^ Y W ^>=1 for (W(l„) WWCO) and (l^.CU 
(3.12) with z Wjbf , . , N ~* being free to be assigned. for I^fO regardless of whether Y^pl or not 

Otherwise* z^*^ " tf{ is preassigned to zero or the Then the recovered feasible solution can be written as 
corresponding arc is not allowed in the feasible set of arcs. 
To ensure that a feasible solution exists, we now only need 
ensure that the variables z /w0 o* - * are free to be assigned 
for l k+L =Q+L . . . JLjfc^ with the corresponding cost coeffi- 
cients being well-defined. (Recall that all variables of the 
form z, 0 0 "-* +1 and % 0 0 . o""^ 1 are free (to be ^ 
assigned) and all coefficients of the form C/ 0 and ^ ^W* and Lower Bound 7116 bound 0D ^ 

Co . . . o^o . . . o"^ 1 are well-defined for 1,=0, .' . X, and y= feasible solutioD is ^ ven b ? 

0. . . . M p for p=k+l, . . . J$.) This is accomplished as V^y^L t .c t t % , " (3.21) 

follows: If the cost coefficient y k y M yo . . . 0 ^ +L is 1 A 1 v 1 " 

well-defined and z, (l v (/ v 0 " lk * +i is free, then define 25 and the lower by ® N f . . . ,u*) for any multiplier value 

cw> . . o"-*=c, ( J Z " . c/^ 1 with z, +i0 . . . 0 ^ hang (u 3 u"). In particular, we have 

free. Otherwise, since all variables of the form z uo 0 ^ 

and zo^o . . . o^ 1 are free to be assigned wim the ■ - * - • ■ . WrfrtS VrfO (3.22) 

corresponding^osts being weH-defined, set c, w o . . 3Q where (u a ^ is mu M p lier value, C 3 . . . r*A 

c 'i(W<» ... o . . . o wher< f ^ term 1S maximizer of <b N f u*), ~ N multidimensional assign- 

omitted if l fc (l*+i)=0 and the second, KJ^i(^i>=Q^* ment problem (3.1) and 1 N is the recovered feasible solution 
Wjt+iH) and i k+l (l k+1 )=0< define c 0 0 =c 0 0 + . defined by (3.20). Finally, we conclude with the observation 

i . * , vt o tu j *^ *a i -a that V iV (z)=V 2 (Y) where Y is the optimal solution of (3. 17) 

The last step k=N-2. The description of the algonthm ^ ^ c * struction of z in ^ons (3.18H3.20). 
endswimk=N+2/^ Reuge of Muldpiiers , Since ^ most computationally 

(3.14) with k-N 2 is expensive part of the algorithm occurs in the nmimization 



0 otherwise. 



of 



40 



Myv starts or multipliers close to the optimal are fundamentally 

subject to £ iat-i-i~. Xjv_i ' irnportant for real-time speed. The purpose of this section is 

J **° to demonstrate that the mult^)lier set obtained at stage k£ 1 

l n - x 45 provide good starting values for those obtained at step k+ 1. 

Z %h-iIn*un^Mn* Theorem 3.xx. Let (u 3 , . . . ,u^) denote a multiplier set 

obtained in the maximum of O^yy 3 . . . .u*). Then this 

$ € (0 1} for all In i i N . multiplier set satisfies the string of inequalities 

5Q *„(^ . . . ,0^^-:(« 4 , ■ . • i^S • - • SO^SOjS VJfr 

Let Y be an optimal or feasible solution to this where ^ fee ^ st ste P » ^ rnaximization of ^multi- 
2^iimensional assignment problem, The recovered feasible P Uers " not changed in me remaking steps. Fuithemore 

solution zT is defined by *> ^ s ^ i, U 2 ^ " ' * T > 
J denote a maximizer of 0^ +1 (u^ 2 u w ). Then, we have 

- (3.18) 

1 if ii = u(^.-.w-i), h = hih-N-il hih-N-i) <t> N ^i^~ k+ \ ... , s 

60 

or if (L^_j. i N ) - (0 t 0) ^(i/'- 3 ) £ <ti < Vjv(i) 

0 otherwise. 

Inequalities depend on solution of 2d assignment prob- 
65 lem. 

Said in a different way, the recovered feasible solution z* is Proof: Suppose we have a value of the multipliers 
then defined by (u** 2 , . . . ,11*) obtained in the maximization of 



5.959.574 



77 



*«-*fj(«* +2 «*) s Minijuize «to-*.i(l"~* +, n<* +! «") 

Subjec. to: ^ 

<t+2 "'JV 



where 
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-continued 

'fc+2 It 



10 



Since trie constraint I, , 

*Af^ i+1 =l for »jt +2 =l 
so that one has 



•Mjt2 ts now imposed, the feasible region is smaller. 



N **P 

zz 

/»=i+2 i,=0 

z 



'jt J Jt+l ' *p-i i p+V~'N | 



p=i+2 J /»=i+2^=0 



15 



20 



These need not be the maximizer; however, we do assume 
that we have solved the above minimization problem opti- 25 
maliy to evaluate 4> w ^ 1 (u* +2 , . . . ,u"). Just as in the 
definition of the earlier recovery problem, let UW^ti)4+i 
(Wi)h t+t Lfct -0 De an enumeration of indices {l A a A+1 ) of the 
optimal solution w 2 of problem (3.4) (or the first 2 indices 
of the solution w Nsk+1 of the relaxed problem (3.2) con- 30 
structed in equation (3.5)) such that w ijt(/t ^ ilk ^-1 for 
0Ati)WwM0,O) and (W^Om0.0) for 
1 A+1 =0 regardless of whether Woo 2 =l or not. 

Then, the final value of the objective function can be 35 
expressed as 

^mW** 2 ^Minimize OV^C^ 1 ;*^ 2 , ...,«*) 

Subject to Z^ 2 . . . ^ %, itlWi+1 >i U2 ^i*** , 0«x) 40 

where 



where the fact that <J>/v-* +1 does not depend on the multiplier 
set u* +2 due the added constraint set. Also, the last equiva- 
lence follows from the fact that < J ) A(=A+1 7*"*' 3 . . . M N ) is the 
relaxed problem (3.xx) (with k replaced by k+1) for the 
recovery problem (3.xx) in step k of the above algorithm. 

Continuing this way, one can compute (u 3 u N ) at the first 

step (k=l), fix them thereafter, and perform no optimization 
at the subsequent steps to obtain 

<V« 3 , . . . . . . . . . ^3(^<J>^VV2) 

where <& 2 Y A 

To explain how to improve this inequality, let 

(u 2+ **"-* +i U "J**H) denote a maximizer of 4>„_ k+l 

(u* +2 , . . . ..u^). Then by the same reasoning one has (3.xx) 
as stated in the theorem. Q.E.D. 
Summary of the Lagrangian Relaxation Problem. (APB) 
Given the multidimensional assignment problem 



(3. xxx> 



z 



't+l'jt+2-'N 



45 



Minimize £ e h ..^zi^ 
i t ..ijf 

Subject to £ = 1 ih = 1 M 2) 

= l and p - 2, .... AT- l) 

Z l t -i N e\Q t l\ fix *U J lt ...^ 



(4.1) 



Jf-t+l 



^=A-t-2lp=0 



If now another constrain set, say the (k+2) rt set, is added to 
this problem, one has 

<«** 2 , ...u H ) = Minimize *jv_i+i(r*~* +l ; «** 2 , s 



Minimize 



z 



z < 



^t»"+iHt+if*+i>"*» Z 



Subject to 



z 



1 



where Cq is arbitrarily defined to be zero and is included 

50 for notatloiftal convenience and where the superscript N on 
both c and z is not an exponent, but denotes the dimension 
of the subscripts and the assignment problem as stated in 
(2.5). In relaxed and recovery problems c Q Q N need not be 
zero! In this problem, we assume that all zero-one variables 

55 z i t . . precisely one nonzero index are free to be 

assigned and that the corresponding cost coefficients are 
well-defined. (This is a valid assumption in the tracking 
environment.) Although not required, these cost coefficients 
with exactly one nonzero index can be translated to zero by 

60 cost shifting without changing the optimal assignment. 
Having explained many of the relaxation features, it is 
now appropriate to present the Lagrangian relaxation 
algorithm, which iteratively constructs a feasible solution to 
the N-dimensional assignment problem (4.1). 

65 Algorithm 4.1 (Lagrangian Relaxation Algorithm). To 
construct a high quality feasible solution, enoted by 2*, of 
the assignment problem (4.1), proceed as follows: 



5.959.574 

79 80 

0. Initialize the multipliers (u** 2 . . . . .u"). e.g.. (u** 2 

U^HO. 0). «i "2 "3 (5.1) 



Minimize 



ZZZ< 



For k-1 N-2, do fafafa 

1. For the Lagrangian relaxed problem (3.2) from the 5 & & 

problem (3.1) by relaxing the last (N-k-1) sets of con- Subject t0 ZI 4 im = i: * 1 Ml > 

strain ts. ' 3 

2. Use a nonsmooth optimization technique to solve ^ s,^ = l, = i, ... . 



Maximize {O w _ 4+ ,hi r €R A V tl for p=k+2, . . . Jf with u/^ being 
fixed} (4.2) 



10 



Jj =0*3*0 

,*+2 



where O AMm . j (u /m ' 2 u") is defined by equation (3.2). The ^ € (0 , n f or ail f vy 3 , 

algorithm following problem (3.9) provides one method for 15 
computing a function value and a subgradient out of the 

subdifferential at („« u"). as required in most To ensure ^ a feas * le sol " tion (5J > ^ wa >' s exists for 

. . . , . n a sparse problem, all variables with exactly one nonzero 

nonsmooth optimization techniques. indfiX vajiables of the f or mz iiw ^. and 2^ for y* 

3. Given an approximate or optimal maximizer of (4.2), 2Q 1 M r and p=l*23 are assumed free to be assigned with 

say (u* +2 ti"), let w 2 assignment problem (3.4) the corresponding cost coefficients being well-defined. This 

corresponding to this rnaximizer of <JVWu* +2 u"). assumption is valid in the tracking environment [**refPoa. 

**refPob]. 

4. Formulate the recovery (N-k)-dimensional problem ^ Lagrangian Relaxed Assignment Problem. The three 
(3.13). modified as discussed at the end of subsection (33) dimensional assignment problem (5.1) has three sets of 
for sparse problems. At this stage. 1 N as defined in (3.15) 25 constraints and one can describe the relaxation by relaxing 

contains the alignment of the indices {i y i^J. anv Q f me three sets of constraints, the description here is 

Formulate the final two-dimensional problem (3.17) and based on relaxing the last set of constraints. A (M 3 +l> 
complete the final recovered solution t N as in (3.20). dimensional multiplier vector (i.e., u 3 eR^ 1 ) associated 

- 1 - - , „ . * * * ~ * • %a with the 3 th constraint set will be denoted by u 3 =(u 0 3 , 

/^"STi f unds ^ et ^ def ^ ed1 ^ Ul * . . . jn M *f with u o 3 s0 being fixed for notational 

(3.2) with k=l, let V^) be the objective function value of ^^ence. (The zero multiplier u 0 3 is used for notational 
the N-dimensional assignment problem in equation (2.5) conV enie D ce.) Now, the three dimensional assignment prob- 
correspondingjo ^feasible solution iT of the constraints in lem (5 ^ is; relaxed t0 a two-dimensional assignment prob- 
(2.5), and let zr~ + (2.5). Then, ^ l em incorporating last set of constraints into the objective 

function via the Lagrangian. Although any constraint set can 
QrAu 3 i^^v^^v^f) (43) be relaxed, we choose the last set of constraints for conve- 

nience. The relaxed problem is 



is the desired inequality. 
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M 1 M 2 W 3 (5.2) 

mXJMWTQ " Mimini2e ^ ^ ^ Z Z Z + 

COMMENTS n =o, 2 *=o 13*0 



«3 



Z Z ' 



1. Step 2 is the computational part of the algorithm. 
Evaluating search procedure requires 99% of the 45 — 
computing time in the algorithm. This part uses two dimen- 
sional assignment algorithms, a search over a large number _ ivtinimize yy + iz? -YV 
of indices, and a nonsmooth optimization algorithm. It is the = fa fa fa W3 + ' 3 ~ fa 1 
second part (the search) that consumes 99% of the compu- 
tational time and this is almost entirely parallelizable. 50 OL , . , , 

J ^ Subject to 2, Z ^vzh =1 * 1+1 = 1 M i 

2. In track maintenance, the warm starts for the initial '2=0*3=0 
multipliers for the first step are available. For the relaxation Q ^ 
procedure, initial multipliers are available in step 2 from the Z a Z ^1^3 * u 1 + 1 = ** 2 
prior step in the algorithm. 55 



ij :c0 13 =0 



3. There are many variations on the above algorithm. For 0ne of ±c major steps in the algorithm is the inaximiza- 

example, one can compute a good solution on the first pass ^ of ^ ft Xums out mat ^ piecewise affine 

(b=l) and not perform the optimization in (2) thereafter. This function of the multiplier vector u 3 , so that the maximization 

yields a great solution. Thus one can continue the optimi- eo of <J> 3 Since many of these algorithms require a function 

zation at the first pass, and immediately compute quality value and a subgradient of 3> 3 y 3 A address this problem in 

feasible solutions to the problem. the next subsection. We note, however, that there are other 

ways to maximize <J> 3 subsection addresses but one such 

Lagrangian Relaxation Algorithm for the 3-Dimensional method. 

Algorithm. (ABP) In this section, we illustrate the relaxation 65 Properties Lagrangian Relaxed Assignment Problem 

algorithm for the three dimensional assignment problem. For a function evaluation of <J> 3 y suboptimal) solution of 

Having discussed the general relaxation scheme, this relaxed! problem (5.2) can be constructed from that of a 
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two-dimensional assignment problem. Then, the nonsmooth 
characteristics of <J> 3 followed by a method for computing a 
subgradient. 

Evaluation of <J> 3 . i 2 ) j3=j3(ii.i 2 ) and a new cost 
function c ; : 2 by: 



82 



w 2 



(5.7) 



cf it2 -c^ for <ii,fe)*lO,0) 

13=0 



Given an index pair (ii,i 2 ), j3 need not be unique, resulting x $ 
in the potential generation of several subgradients (5.11). 
(We further discuss this issue at the end of Subsection 5.23.) 
Now, define 



<K« 3 ) = Minimize ^(z 2 ;* 3 ) = £ X^i'2^i'2 " £ 

*! =0*2=0 O =0 



(5.4) 



20 



«2 



Subject to £ ^ j/2 - 1, 1, 



i 2 =o 



25 



Woo 2 =0 if (i lf i 2 H00) and c^+u^ for aU i 3 
Woo 2 =l if (i^MOO) and c^+u^^O for some i 3 . 

Then w 2 is a feasible solution of the problem (5.4) and 
<J> 3 (w 3 ;u 3 )!=<J> 3 (w 2 ;u 3 ). If, in addition, w 3 is optimal for 
(5.2), then w 2 is an optimal solution of (5.4), <J> 3 (w 3 ;u 3 )= 
<t> 3 (w 2 ;u 3 ) and cj> 3 ( u 3 )=0 3 (u 3 ). 

The Nonsmooth Optimization Problem 
An Algorithm for Computing <J> 3 and a Subgradient 
Given u 3 the problem of computing 0 3 y 3 A subgradient of 
<J> 3 y 3 A 

1. Starting with problem (5.1), form the relaxed problem 
(5.2) 

2. To solve (5.2) optimally, define the two dimensional 
assignment problem (5.4) via the transformation (5.3) 

3. Solve the two-dimensional problem (5.4) optimally. 

4. Reconstruct the optimal solution, say w 3 e 7 A equation 
(5.6) as in Theorem 5.1. 

5. Then 



Z^ = l > 1 

sf lf2 e{0,l\ for all i u h- 



M x M 2 M3 
i^=0 <2=0 (3=0 



(5.10) 
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As an aside, two observations are in order. The first is that 
the search procedure needed for the computation of the 
relaxed cost coefficients in (5.3) is the most computationally 
intensive part of the entire relaxation algorithm. The second 35 
is that a feasible solution z 3 of a sparse problem (5.1) yields 
a feasible solution z 2 of (5.4) via the construction 



6. A subgradient is given by substituting w 3 objective 
function (5.2), erasing the minimum sign, and taking the 
partial with respect to u^ 3 . The result is 



1, if = 1 for ('l* l2 ' * 3); 



40 



0, otherwise. 



4 



<5mf 3 



*l=0i2=0 



(5.U) 



for k = 1, 



, Af>. 



Thus the two-dimensional assignment problem (5.4) has 
feasible solutions other than those with exactly one nonzero 
index if the original problem (5.1) does. 

The following Theorem 5.1 states that an optimal solution 
of (5.2) can be computed from that of (5.4). Furthermore, if 
the solution of either of these two problems is e e then so is 
the other. The converse is contained in Theorem 5.2. 

Theorem 5.1. Let (0 2 dimensional assignment problem 
(5.4) and define w 3 by 

w^ 3 =0 if i 3 *j 3 and (i„i 2 )*(0,0); 

Woo^lifcoo^-Hi^. (5.5) 

Then w 3 is a feasible solution of the Lagrangian relaxed 
problem and <1> f 3 A O f 3 A e e 2 optimal for the two 
dimensional problem, then w 3 is an optimal solution of the 
relaxed problem and <E> 3 (u 3 )=<I> 3 (u 3 ). 

With the exception of one equality being converted to an 
inequality, the following theorem is a converse of Theorem 
5.1. 

Theorem 5.2. Let w 2 be a feasible solution to problem 
(5.2) and define w 2 by 



The Recovery Procedure. The next objective is to explain 

45 a recovery procedure, i.e., given a feasible (optimal or 
suboptimal) solution w 2 of (5.4) (or w 3 of (5.2) constructed 
in Theorem 5.1). generate a feasible solution z 3 of equation 
(5.1) which is close to w 2 in a sense to be specified. We first 
assume that no variables in (5. 1) are preassigned to zero; 

50 this assumption will be removed shortly. The difficulty with 
the solution w 3 in Theorem 5.1 is that it need not satisfy the 
third set of constraints in (5. 1). (Note, however, that if w 2 is 
an optimal solution for (5.4) and w 3 , as constructed in 
Theorem 5.1, satisfies the relaxed constraints, then w 3 is 

55 optimal for (5.1).) The recovery procedure described here is 
designed to preserve the 0-1 character of the solution w 2 of 
(5.4) as far as possible: If w Iifz 2 =l and i t *Q or i 2 *Q* the 
corresponding feasible solution z 3 of (5.1) is constructed so 
that z lil2 , 3 3 = 1 for some (i 1 .i 2 i 3 ) . By this reasoning, variables 

60 of the form z^ 3 , can be assigned a value of one in the 
recovery problem only if Woo 2 =l. However, variables Zoo, 3 3 
will be treated differently in the recovery procedure in that 
they can be assigned either zero or one independent of the 
value Wqq 2 . Th* s increases the feasible set of the recovery 

65 problem, leading to a potentially better solution. 

{(iiCy^*?.)}^-^ be an enumeration of indices 
Vh>h} °f w 2 ( or me fas* 2 indices of w 3 constructed in 
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equation (5)) such that w (lW ft) 2 =l for (i^UM^l^MCXO) Let Y be an optimal or feasible solution to this 

and (ij(l 2 )i 2 (l 2 ))==(0.0) for 1 2 =0 regardless* of whether 2-dimensional assignment problem. The recovered feasible 

Woo 2 =l or not. (The latter can only improve the quality of the solution z 3 is defined by 
feasible solution.) 

Next to define the two-dimensional assignment problem 5 ^ J i if h « UikX h = k(k) and y^ = l; (515) 

that restores feasibility, let ^ = 1 0 otheiw ise. 

fa ^ 0 -"^ (5 - 12) 

Said in another way. let { (l 2 (l 3 )i 3 (l 3 ) }t^=0 be an enumera- 
Then the two-dimensional assignment problem that restores 10 tion of indices {l 2 i 3 } of Y such mat Y^^^l for (U(h)* 
feasibility is i 3 (l 3 M0.0) and ~(1 2 (1 3 ) j 3 (1 3 M0.0) for 1 3 A regardles's of 

whether Y C)0 =1 or not. Then the recovered feasible solution 
h **a (5.14) can be written as 



/ 2 =Oi 3 =o 15 
Subject to £ 4 3 = U 2 = 1 I 0 othfi ™se 



1 if ii = iiUnl h = i 2 (/i2\ i 3 = hih) {5A6) 



i 3 =o 



£4,3-1,13-1 Mj, 



20 The upper bound on the feasible solution is given by 



/ 2 «0 



Mi jw 2 W3 (5.17) 
€ {0, 1} for all t ly i 5 . V 3 U 3 ] - £ Z Z 4*3* 

25 

The next objective is to remove the assumption that all ^ ^ lowef by <j> Y 3 A a A 
cost coefficients are defined and all zero-one variables are ^ particular we have 
free to be assigned. We first note that the above construction " 

of a reduced dimension assignment problem (5.13) is valid &3(i?)£& 3 {i?y£v${£)£V 3 {?) 

as long as all cost coefficients c 2 are defined and all zero-one 30 where u 3 fc multiplier value< $ A 3 0 y 3 Ae ? assign _ 
vanables in z are free to be assigned. Modifications are ment blem and ? is ^ recovered feasible solution 
necessary for sparse problems. If the cost coefficient defined by (5 16) 

*! well-defined and the zero-one variable 1^.^ for ^ Multidimensional Assignment 

z ^h->h 15 fr f 10 ^^''ff t0 zero Qr 3°/ ,e - *?» defi f Problem.(ABP) In this section, we briefly describe other 
c ^"- c 'i«2teft> 5 m ( 5 - 1 2-> wrthz . 1 <« i )««i, being free to be 35 possible relaxations and their impUcations. These include 
assigned. Otherwise. z h is preassigned to zero or the ^ programming relaxations and the corresponding 
corresponding arc is not allowed in the feasible set of arcs. duals 

To ensure that a feasible solution exists, we now only need Recall ^ one stam with &c definition of ^ zero . one 
ensure that the variables z^ 0 are free to be assigned for yariabie j and CQSt coefficients t0 form me problem 

1 2 =0,1 with the corresponding cost coefficients being 40 @ 

weU-defmed. (RecaU that all variables of the form z li00 3 . wher£ ig ^ defined to ^ ^ Here , each 

z 3 and z^ 3 are free (to be assigned) and all coefficients ^ rf sums in ^ constraints repr esents the fact that each 
of the corresponding form c 4i00 and c^ .) This is non _ dummy report occurs exactly once in a "track of data", 
accwjp^asf^ is Qne can modify ^ formulation to include multiassign- 

weU-defined and z^ hl0 is free, then define c^ 45 men ts of one, some, or all the actual reports. The assignment 

m with Zjb0 being free. Otherwise since all vanables of blem ^ i$ ch d accordmgly . For exara ple, if z * is 

the form i* and z 3 are free to be assigned with the tQ ^ ^ d no more ^ ^ or n0 less ^ * 
corresr^nding costs being well-defined, set c»^<uoo + ^ ^ ^ ti=r in ^ constraint (2 5) is A d £ 
co^p 3 where the first term is omitted if 1 (1 2 >=0 and the «^ _ > respectively. Modifications for group tracking 
second if i 2 (l 2 )=0. For i^M) and i 2 (l 2 )=0. define ^=5 0and mu itiresolution features of the surveillance region will 
c °°° ' be addressed in future work. In making these changes, one 

The final recovery problem. The recovery problem for the must pa y careful attention to the independence assumptions 
case N=3 is that need not be valid in many applications. 

Again, the recent work of Poore and Drummond signifi- 
.. , 1 ; (5.14) 55 cantly extends the formulation of the track maintenance and 

2j ^3 initiation to new approaches. The discussions of this section 

apply equally to those formulations. 
Vz? Ln The Linear Programming Relaxation. In the linear pro- 

u ^ 2j^i^ " * 2 ~ *"-*^ 2 ' CTammins relaxation, one reolaces the zero-one variables 
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gramming relaxation, one replaces the zero-one variables 
60 constraints 
h 

with the constraint 

4a € 1} for /2 ' i3 ' 65 OSz^ . . . far «n . . . v (6.2) 

Then, the problem (6.1) can be formulated as a linear 
programmiing problem with the constraint (6.2) in (6.1) 
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replaced by (6.3) with a special block structure to which the 
Dantzing-Wolfe decomposition is applicable. Of course, 
after solving this problem, one must now recover the zero- 
one character of the variables in (6.1) and there are many 
ways to do this, such as using the two dimensional assign- 
ment problems. Commercial software is also available. 

The Linear Programming Dual and Partial Lagrangian 
Relaxations. Given the linear programming relaxation, one 
can formulate the dual problem or the partial Lagrangian 
relaxation duals with respect to any number of constraints. 
In particular, this is precisely what is done in Section 3 on 
the Lagrangian relaxation algorithm presented. The much 
broader class of algorithms provided in the U.S. Pat. No. 
5,537,119 (Aubrey B. Poore, Jr., Method and System for 
Tracking Multiple Regional Objects by MultiDimensional 
Relaxation, U.S. Pat. No. 5,537,119, filed Mar. 14. 1995, 
issue date of Jul. 16, 1996) can be modified to remove the 
zero-one character when one relaxes M sets of constraints to 
an N-M dimensional problem and recovers with an N-M-s-1 
dimensional problem. This avoids the problems associated 
with the NP-hard N-M and N-M+l dimensional problems. 
However, to restore the zero-one character, one can do it 
sequentially with an assignment problem or with one of the 
many zero-one rounding techniques. These formulations are 
easy to work out and thus the details are omitted. 

The complete dual problem is another way of solving the 
LP problem and may indeed be more efficient in certain 
applications. In addition, the solution of this dual problem 
may provide excellent initial approximations to the multi- 
pliers for Lagrangian relaxations. 

Multisensor Resolutions Problems and Their Solution via 
Linear Nonlinear Assignment Problems and Inequality Con- 
straints. (ABP) The mathematical development for inequal- 
ity constraints is understood by Aubrey B. Poore, but this 
work has not been written up. This feature is important in 
combining information from sensors with different resolu- 
tions such an IR and Radar. 

Formulation and Algorithms for the Distributed Sensor 
Tracking: Track-to-Track Correlations. (ABP) 

Hot Starts for Track Maintenance and Initiation: Bundle 
Modifications (ABP) 

Thus reuse of multipliers and the first proof that this reuse 
is actually valid was presented in this section. The reuse in 
track maintenance is demonstrated and discussed in the 
work of Poore and Druramond [xx]. The only item left is the 
modification of the bundle of subgradients for the use with 
the multiplier values as one goes through the recovery 
problem as in Section 3.6 or as one moves the window in 
track maintenance. It is this aspect of the nonsmooth opti- 
mization that adds an order of magnitude to the speed of the 
algorithms. This work is based on both a mathematical proof 
as in Section 3.6 as well as computational experience and 
heuristics. 

Adaptive Window Size Selection (ABP) 

These data association algorithms, based on the multidi- 
mensional assignment problem, range from sequential pro- 
cessing to many frames of data processed all at once! The 
data association problem for multisensor multitarget track- 
ing is formulated using a window sliding over the frames of 
data from the sensors. This discussion focuses on the work 
of Poore and Drummond and not on the formulations in 
Section 2. Firm data association decisions for existing track 
are made at, say frame k, with the most recent frame being 
k+M. Decisions after frame k are soft decisions. Reports not 
in the confirmed tracks are used to initiate tracks over frames 
numbered k-I to k+M. 
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The length of these windows varies from sequential 
processing, which corresponds precisely to 1=0 and M=l. to 
user defined large values of I and M. (In the case of 
sequential processing, we also have a temporary track file of 

5 dynamically feasible tracks, but incorrect data association.) 
There is a marked change in performance over this range. As 
the two windows become exceedingly long, there is little 
statistically significant information gained and indeed per- 
formance can degenerate slightly. Thus, one must manually 

10 find the correct window lengths for performance in a given 
scenario and the algorithms do not change for a changing 
environment. Thus we propose an adaptive method for 
adjusting the window lengths. (The method has been highly 
successful for selecting the correct window length for mul- 

15 tistep methods in ordinary differential equations.) 

1. Compute the solution for different window lengths that 
overlap and differ by one or two frames. 

2. Compare the solution quality, i.e., the value of the 
;o objective function, for two adjacent window lengths. 

3. If there is a predefined improvement in either direction, 
we then, for stability, repeat the exercise for a shorter or 
longer on the first try. If there is consistent information, we 
adjust the window size in the indicated direction. This can be 

l5 given a well defined mathematical formula in terms of the 
assignment problems of different dimensions. 
Algorithm Enhancements due to Data Structures (ABP) 
Construction and enumeration of the Best K- Solutions. 
New Nonsmooth Optimization Techniques. 

30 Sensitivity Analysis 

For sequential processing in which the two dimensional 
assignment algorithm is used, one can use the LP sensitivity 
analysis theory and easily obtain the corresponding answers. 

3 5 For the multiframe processing, the optimal solution is not 
available; however, there are still two approaches one can 
use. First, the basic algorithm can perform the same sensi- 
tivity analysis at each stage (loop) of the algorithm as is done 
in the two dimensional assignment problem since the evaiu- 

40 ation of function & N - k +i is equivalent to a two dimensional 
assignment problem Alternately, one could use an LP relax- 
ation of the assignment problem and base the sensitivity on 
the resulting LP problem. We currently see this as an 
important step in finding even better solutions to the assign- 

45 ment problem if so desired. 

New Auction Algorithms (ABP & PS). In this section we 
present a new auction algorithm for the two dimensional 
assignment problem. 
An important step in solving N-dimen'sional assignment 

50 problems for N^3, is finding the optimal solution of the 
2-dimensional assignment problem. In particular, we wish to 
solve the 2-dimensional problem which includes the zero- 
index. T3,picaUy this problem can be thought of as trying to 
find either a minimum cost or maximum utility in assigning 

55 person to objects, tenants to apartments, or teachers to 
classrooms. We will follow the work of Bertsekas and call 
the first index set persons and the second index set objects. 
The statement of this problem is given below (1) when the 
number of persons is m and the number of objects is n. 

60 

» ■ (1) 

Minimize ^Z/'j^ 

N 

65 Subject to ^Xij = 1 for all i = 1 T ... , m 



87 

-continued 



5.959.574 



88 



cyHif+Uj 2 ^ for all (iJ)£S, 



There are a couple of assumptions which we will make 
about (1). First of all, we assume that c l0 and c 0j are 
well-defined and the corresponding variables % l0 Xo y free to 
be assigned Second if a cost c iy happens to be undefined, 
then the corresponding variable % tJ % to 0. In effect if c v is 
undefined, we simply remove this cost and variable from 
inclusion in the problem, 

Notice that there are no constraints on the number of times 
person 0 and object 0 can be assigned. But notice that the 
first constraint requires that each non-zero person i must be 
assigned exactly once. Similarly, the second constraint 
forces each non-zero object j to be assigned exactly one 
time. Finally, the last constraint gives an integer solution, 
although we will see shortly that this constraint can be 
relaxed to admit any solution % tJ ^ 0. One reason for requir- 
ing that all of the costs of the form c (0 and c 0y be defined is 
so that we are guaranteed a feasible solution exists for the 
given problem. 

Relaxation of One Constraint 

Relaxing the last constraint, we obtain the following 
problem: 



10 



m n \ n " \ t m 

,=0 U=0 



£ Zi>*- + ^-z*5 



j=0 jfQ 



;=0 



Subject to 2^ = 1 for all / = 1, ... , m 
XU € {0, 1} 



For i = 0, xoj-- 



1 if coj + t^ <0 



[0 otherwise 

For i * 0, 
All of fs are assigned only one time. 



{1 if y-argmin{c# + tfj |* = 0, l t ...,«} ] 
0 otherwise j 



The relaxation of the first constraint is analogous and would 
lead to similar results with i and j exchanged and the 
introduction of the multiplier u 1 . 

Relaxation of Both Constraints 

This time we will relax both constraints and using duality 
theory obtain a relationship between the multipliers u 1 and 
u 2 . 

a bunch of equations, theorems and the such here. 

From the above relaxation, we arrive at the following 
complementary slackness conditions which have been 
relaxed by a small value e. 

Definition: An assignment S and a multiplier set (u\u 2 ) are 
said to satisfy e-complementary slackness (e-CS) if 



c^+uf+ufZ-e for all (i,j>LA. 

Forwar d Auction Algorithm 

(1) Select any unassigned person i 

(2) Determine the following quantities: 

jf=arg min \ c & +u k 2 \teA([)} 
KV=min \c A +u k 2 UxA(t), k*j,} 



In the selection of j, above, if a tie occurs between 0 and 
is any non-zero index k. then select j, as k Otherwise, if there 
is a tie between two or more non-zero indices, the choice of 
j, is arbitrary. Also if A(i) consists of only one element, then 
set w l -=oo. 

(3) Update the multipliers and the assignment: 
If jr=0, then 

(a) Add (i,0) to the assignment. 

(b) Update u/: -f-c l0 . 
jf*0, then 

(a) Add (Lj,) to the assignment 

(b) Remove (i'.j z ) from the assignment if j, was previously 
assigned. 

(c) Update u,. 2 : =u y 2 4<w -v,)+€=w -c y +e. 

(d) Update u, 1 : ^-(c^ii^)--^^ 
Reverse Auction Algorithm 

(1) Select any unassigned object j, such that c^+u/zsO. 

(2) Determine the following quantities; 
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ipar.g min {c^+u^iheBif)} 
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This lower bound is achieved and the second constraint in 
the original problem is satisfied if the following conditions 
are met: 
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fj=mm \c kj +u k l \k*B(j), k*i,} 

In the selection of i^ above, if a tie occurs between 0 and 
any non-zero index L then select j, as k. Otherwise, if there 
is a tie between two or more non-zero indices, the choice of 
j f is arbitrary. Also if B(j) consists of only one element, then 
set y y oo. 

(3) Update the multipliers and the assignment: 
If i=0. then 

(a) Add (0,j) to the assignment. 

(b) Update u/: =-Co r 
If i/O. then 

(a) Add (i r j) to the assignment 

(b) Remove (i r j') from the assignment if (i y . j) was previ- 
ously assigned. 

(c) Update u^ 1 : ^^y^p^^^Cy+e 

(d) Update u/: ^c+u^-yj-e. 
Combined Forward/Reverse Auction Algorithm 

1. Assume that u 2 is given as an arbitrary multiplier. 

2. Adjust the value of u 2 for each object j as follows: 

If c Oy -Hi/<0, then set u/: =~c 0J . 



3. Run iterations of the Forward Auction Algorithm until 
all persons become assigned. 

4. Run iterations of the Reverse Auction Algorithm until 
65 all of the objects become assigned. 

Note at the completion of the Forward auction step we 
have the following conditions satisfied: 
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c,y+u/+u 2 =0 for all (i.j)eS. 

c^-m/Sir"^"* 2 } for all (ij)eS. 

Thus we can prove the following proposition. 
Proposition: If we assume that c o/ Hi/^0 at the start of the 
Forward Auction Algorithm and all of the persons are 
assigned via a forward step, then we have: 



for all (i j)eA. 
Ctf+uj+uj 2 ^) for all (ij)eS. 



Optimality of the Algorithm 

Theorem: e-CS preserved during every forward and 
reverse iteration. 

Theorem: If a feasible solution exists, then the resulting 
solution is with me of being optimal for the Combined 
Forward Reverse Algorithm, 

Implementation Specifics 

Parallelization 

Here are but a few comments. 

Although the algorithm appears to be serial in nature, its 
primary computational requirements are almost entirely 
parallelizable. Thus parallelization is planned. 

Step 2 is the computational part of the algorithm. Evalu- 
ating <&N-k+i XXX procedure requires 99% of the computing 
time in the algorithm. This part uses two dimensional 
assignment algorithms, a search over a large number of 
indices, and a nonsmooth optimization algorithm. It is the 
second part (the search) that consumes 99% of the compu- 
tational time and this is almost entirely parallelizable. 
Indeed, there are two dimensional assignment solvers that 
are highly parallelizable. Thus, we need but parallelize the 
nonsmooth optimization solver to have a reasonably com- 
plete parallelization. 

If a sensitivity analysis is desired or if one is interested in 
computing several near-optimal solutions, a parallel proces- 
sor with a few powerful processors and good communica- 
tion such as on the Intel Paragon would be most beneficial. 

The foregoing discussion of the invention has been pre- 
sented for purposes of illustration and description. Further, 
the description is not intended to limit the invention to the 
form disclosed herein. Consequently, variation and modifi- 
cation commensurate with the above teachings, within the 
skill and knowledge of the relevant art, are within the scope 
of the present invention. The embodiment described here- 
inabove is further intended to explain the best mode pres- 
ently known of practicing the invention and to enable others 
skilled in the art to utilize the invention as such, or in other 
embodiments, and with the various modifications required 
by their particular application or use of the invention. It is 
intended that the appended claims be construed to include 
alternative embodiments of the invention to the extent 
permitted by the prior art 
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What is claimed is: 

1. A melhod for tracking a plurality of objects, compris- 
ing: 

repeatedly scanning a region containing a set consisting of 
one or more moving objects and generating N sequen- 
tial images or data sets of said region, a plurality of 
observations in said images or data sets providing 
positional information for objects in said set; 

determining a plurality of tracks, at least one track for 
each object in said set; 

determining a plurality of costs, wherein each cost is for 
assigning one of said observations to one of said tracks; 

defining a linear programming problem: 
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Minimize 



Z<v 



Subject To *y Zh 
..i N 

Z *v 



■iff 



= 1 ('I = 1, 



. Mi) 



• V = 1 ('2 = 1 Af 2 ) 



= 1 



(i p = 1 M p and p = X I) 



Z 



0<Zi v „ lN < 1 for all h. 



wherein each c (i tff is included in said plurality of costs, 
each M,, i=l, \ . . *,N, being one of: (a) a number of 
observations in an I th image or data set of said N sequential 
images or data sets; (b) a sum of a number of tracks in said 
plurality of tracks, and a number of said observations in the 
i"' image or data set not assigned to one of said tracks; and 
(c) a number of tracks in said plurality of tracks; 

solving said linear programming problem for values of 

z . , for each il . . . iN; 
determining a value z (i . . . ifl in {0,1} for each il . . . iN 
corresponding to each z t; i^, wherein said values 
zil . . . iN provide an optimal or near optimal solution 
to said linear programming problem; 
taking one or more of the following actions based on said 
optimal or near-optimal assignment of said plurality of 
points to said plurality of tracks: 
sending a warning to aircraft or a ground or sea facility, 
controlling air traffic, 

controlling anti-aircraft or anti-missile equipment, 
taking evasive action, 

working on one of said one or more objects, surveilling 
one of said one or more objects. 



